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Abstract 

In this thesis we study graph algorithms, both in sequential and parallel contexts. In 
the following outline of the thesis, algorithm complexities are stated in terms of the number 
of vertices n, the number of edges m, the largest absolute value of capacities U, and the 
largest absolute value of costs C. 

In Chapter 1 we introduce a new approach to the maximum flow problem that leads 
to better algorithms for the problem. These algorithms include an O(nm\og(n 2 /m)) time 
sequential algorithm, an 0(n 2 log n) time parallel algorithm that uses 0(n) processors and 
0(m) memory, and both synchronous and asynchronous distributed algorithms. 

Chapter 2 is devoted to the minimum cost flow problem, which is a generalization 
of the maximum flow problem. We introduce a framework that allows the generalization 
of the maximum flow techniques to the minimum-cost flow problem. This framework al- 
lows us to design efficient algorithms for the minimum-cost flow problem. We exhibit 
0(nm log(n)log(nC)), 0(n 5 / 3 m 2 / 3 log(nC)), and 0(n 3 log(nC)) time sequential algorithms 
as well as parallel and distributed algorithms. 

In Chapter 3 we address implementation of parallel algorithms through a case-study of 
an implementation of a parallel maximum flow algorithm. Parallel prefix operations play 
an important role in our implementation. We present experimental results achieved by the 
implementation. 

Parallel symmetry-breaking techniques are the main topic of Chapter 4. We give an 
0(\g*n) algorithm for 3-coloring a rooted tree. This algorithm is used to improve several 
parallel algorithms, including algorithms for A+1-coloring and finding maximal independent 
set in constant-degree graphs, 5-coloring planar graphs, and finding a maximal matching 
in planar graphs. We also prove lower bounds on the parallel complexity of the maximal 
independent set problem and the problem of 2-coloring a rooted tree. 
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Introduction 



Advances in computing technology have enabled computers to solve a wider class 
of problems. One way to extend this class of problems is to build more powerful 
computers. Another way is to design better algorithms. Both approaches have their 
advantages and disadvantages. Fast hardware is essential for many applications, 
such as these that require real-time computations. Furthermore, any program runs 
faster on a faster machine. On the other hand, if an algorithm takes too much 
time to solve real-life problems due to high (e.g., exponential) running time, the 
combinatorial explosion will prevent the algorithm from solving these problems 
even if computers become thousands of times faster. 

Parallel computers are one of the most promising recent developments in high- 
performance hardware. Parallel machines with tens of thousands processors are 
already available commercially, and machines with millions of processors can be 
built with today's technology. The class of problems for which efficient parallel 
algorithms are known, however, is much smaller then the class of problems for 
which efficient sequential algorithms are known. The implementation of parallel 
algorithms is also not as well understood. 

This thesis addresses the issue of designing efficient sequential and parallel al- 
gorithms for graph-theoretic problems. These problems are important because of 
their applications, both outside of and in the field of computer science. For this 
reason, graph algorithms have become one of the most studied fields of theoreti- 
cal computer science. We shall see, however, that even algorithms for the classical 
problems in the area can be significantly improved. 

This thesis has four chapters. In Chapter 1 we study the maximum flow problem, 
a classical problem in the theory of network flows. We introduce a new approach 
to the problem and use this this approach to design better sequential, parallel, and 
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distributed algorithms for the problem. In particular, we exhibit 0{nm\og(n 2 / m)) 
sequential algorithm for the problem, improving the previous upper bound. A 
parallel version of our algorithm runs in O(n 2 logn) time using 0(n) processors 
and 0(m) memory. These time and processor bounds are as good as for those the 
best of the previously known algorithms, and the memory bound is better, making 
our algorithm much more practical. We also exhibit and analyze synchronous and 
asynchronous distributed algorithms. 

In Chapter 2 we study another classical problem from the theory of network 
flows: the minimum-cost flow problem. We show how maximum flow techniques (in- 
cluding the techniques introduced in Chapter 1) can be extended to the minimum- 
cost flow problem. We introduce a new framework for the minimum-cost flow prob- 
lem and exhibit algorithms that significantly improve the known complexity bounds. 
Our sequential algorithms achieve 0(min(nm log re, n 5 / 3 m 2 / 3 , n 3 ) log(nC)) running 
time. A parallel version of our algorithm achieves 0(n 2 log(n)log(nC)) running 
time using 0{n) processors and 0(n 2 ) memory. A different parallel algorithm uses 
only 0(m) memory and 0{n) processors. The parallel running time bound we can 
prove for this algorithm is 0(n 3 log(n)log(nC)); however, we conjecture that the 
running time is in fact 0(n 2 log(n)log(nC)). 

How does one implement parallel algorithms? Are parallel algorithms really 
faster then the sequential algorithms? These questions are addressed in Chapter 
3, where we describe a parallel implementation of our maximum flow algorithm 
that uses parallel prefix operations as primitives. In this chapter we also present 
experimental results achieved by the implementation. 

Our study of design and implementation of parallel network flow algorithms 
motivates the search for general techniques for the design of parallel algorithms. 
Chapter 4 introduces efficient symmetry-breaking techniques which are very impor- 
tant in the context of parallel computation. We show how to apply these techniques 
to design better parallel algorithms. Our main result is an 0(lg*n) parallel time 
algorithm to 3-color a rooted tree. This algorithm is used to improve parallel algo- 
rithms for several graph-theoretic problems, including the problem of 5-coloring a 
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planar graph, the problem of (A + l)-coloring and finding a maximal independent 
set in a constant-degree graph, and a few other problems. All these algorithms use 
a linear number of processors. We also show that the 2- coloring of a rooted tree 
requires fi(logra/ log log rc) time if a polynomial number of processors is used, and 
prove the same lower bound for the maximal independent set problem on general 
graphs. 






fees?' 



si 



*»•: 






m 









\ ' 



. : i^'^ : --r-.^^ 


; : - 






?'.;• • 


r 


W" 


,-_ 










^*\ 




- ##>- '. 








' .<$L.. . 






'.Jftjjfe! 



Chapter 1 



The Maximum Flow Problem 



1.1 Introduction 



In this chapter we introduce a new approach to the maximum flow problem and use 
this approach to design better sequential and parallel algorithms for the problem. 

The problem of finding a maximum flow in a directed graph with edge capacities 
arises in many settings in operations research and other fields, and efficient algo- 
rithms for the problem have received a great deal of attention. Extensive discussion 
of the problem and its applications can be found in the books of Ford and Fulker- 
son [19], Even [16], Lawler [52], Papadimitriou and Steiglitz [60], and Tarjan [71]. 
Table 1.1 summarizes polynomial-time algorithms for the problem. Time bounds 
are stated in terms of the number n of vertices, the number m of edges, and in two 
cases in terms of an upper bound U on the edge capacities (assumed in these cases 
to be integers). 

The first maximum flow algorithm, due to Ford and Fulkerson [20], works by 
finding augmenting paths. The algorithms in the table are variations of the Ford- 
Fulkerson algorithm, incorporating the observation of Edmonds and Karp [15] that 
augmenting along shortest paths leads to a polynomial-time algorithm (algorithm 
1). To further improve the efficiency, Dinic [14] proposed a method to find all 
shortest augmenting paths in one phase. Algorithms 2-11 use Dinic's method. Al- 
ii 
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CHAPTER 1 . THE MAXIMUM FLOW PROBLEM 



# 


Date 


Discoverer 


Running Time 


References 


1 


1969 


Edmonds and Karp 


0(nm 2 )) 


[15] 


2 


1970 


Dinic 


0(n 2 m)) 


[14] 


3 


1974 


Karzanov 


0(n 3 ) 


[47] 


4 


1977 


Cherkasky 


0(n 2 m}l 2 ) 


[12] 


5 


1978 


Malhotra, Pramodh Kumar, 
and Maheshwari 


0(n 3 ) 


[56] 


6 


1978 


Galil 


0{n 5 l 3 m 2 l 3 ) 


[27] 


7 


1978 


Galil and Naamad; 
Shiloach 


0(nm(logn) 2 ) 


[29] 

[64] 


8 


1980 


Sleator and Tarjan 


0(nm log n) 


[67,68] 


9 


1982 


Shiloach and Vishkin 


0{n 3 ) 


[66] 


10 


1983 


Gabow 


O(nmlogU) 


[26] 


11 


1984 


Tarjan 


0(n 3 ) 


[72] 


12 


1985 


Goldberg 


0(n 3 ) 


[33] 


13 


1986 


Goldberg and Tarjan 


0(nm log(n 2 /m)) 


[37] 


14 


1986 


Ahuja and Orlin 


0(nm + n 2 \ogU) 


[1] 



Table 1.1: Polynomial-time algorithms for the maximum flow problem. Algorithm 13 
is presented in this chapter. 



gorithms 12-14 are based on the approach described in this paper. 

There is no clear winner among the algorithms in the table that are based 
on Dinic's method. Algorithms 3, 5, 9 and 11 are designed to be fast on dense 
graphs, and algorithms 4, 6, 7, 8, and 10 are designed to be fast on sparse graphs. 
For dense graphs, the best known bound of 0(n 3 ) was first obtained by Karzanov 
[47]; Malhotra, Pramodh Kumar, and Maheshwari [56] and Tarjan [72] have given 
simpler 0(n 3 )-time algorithms. For sparse graphs, Sleator and Tarjan's bound of 
0(nm log n) [67,68] is the best to date. For a small range of densities (m between 
n(n 2 /(logn) 3 ) and 0(n 2 )), Galil's bound of 0(n 5 l 3 m 2 ' 3 ) [27] is best. For sparse 
graphs with integer edge capacities of moderate size, Gabow's scaling algorithm 
[26] is best. Among the algorithms based on Dinic's method, the only parallel 
algorithm is that of Shiloach and Vishkin [66]. This algorithm has a parallel running 
time of 0(n 2 logn) but requires 0(nm) space. Vishkin (private communication) has 
improved the space bound to 0(n 2 ). Our work has been motivated by the Shiloach- 
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Vishkin algorithm. 

In this chapter we present a different approach to the maximum flow problem, 
which is the basis for algorithms 12-14 in the table. Our method uses Karzanov's 
idea of a preflow. A preflow is like a flow except that the total amount flowing into 
a vertex can exceed the total amount flowing out. During each phase, Karzanov's 
algorithm maintains a preflow in an acyclic network. The algorithm pushes flow 
through the network to find a blocking flow, which determines the acyclic network 
for the next phase. Our algorithm abandons the idea of finding a flow in each phase, 
and also abandons the idea of global phases. Instead, our algorithm maintains a 
preflow in the original network and pushes local flow excess toward the sink along 
what it estimates to be shortest paths in the residual graph. This pushing of flow 
changes the residual graph and paths to the sink may become saturated. Excess 
that cannot be moved to the sink is returned to the source, also along estimated 
shortest paths. Only when the algorithm terminates does the preflow become a 
flow, and then it is a maximum flow. 

The algorithm is simple and intuitive. It has natural implementations in se- 
quential and parallel models of computation. We present a simple sequential imple- 
mentation that runs in 0(n 3 ) time and a more complicated sequential implemen- 
tation that uses the dynamic tree data structure of Sleator and Tarjan [68,69,71] 
and runs in 0(nm\og(n 2 / m)) time. The latter bound matches the best known 
bounds as a function of n and m for both sparse and dense graphs and is better 
than known bounds on graphs of intermediate density. We present a parallel ver- 
sion of the algorithm running in 0(n 2 logn) time using 0(1) words of storage per 
edge. This matches the time bound of the Shiloach- Vishkin algorithm, but our im- 
proved space bound allows implementation on a model of distributed computation 
in which the amount of space per processor at a vertex is bounded by the vertex 
degree. Recently, Ahuja and Orlin [1] used the approach described in this chapter to 
develop an 0(nm + n 2 log U) algorithm for the problem, improving Gabow's bound 
of 0{nm log U) [26]. 

This chapter contains six sections in addition to the introduction. Section 1.2 
describes a generic version of the algorithm. Section 1.3 proves its termination and 
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correctness. Section 1.4 refines the algorithm to produce an 0(n 3 )-time sequential 
implementation. Section 1.5 adds the use of dynamic trees and thereby improves 
the sequential time bound to O(nmlog(n 2 /m)). Section 1.6 discusses efficient dis- 
tributed and parallel implementations. Section 1.7 contains some concluding re- 
marks and open problems. 

The approach presented in this section has been pioneered by the author [33]. 
The version presented here generalizes and improves the original results. These 
improvements represent joint work with Tarjan [37]. 

1.2 A Generic Maximum Flow Algorithm 

Let G = (V,E) be a directed graph with vertex set V and edge set E. We shall 
denote the size of V by n and the size of £ by m. For ease in stating time bounds 
we assume m > n — 1 > 4. For a pair of vertices v and w we define the distance 
d G (v, w) from v to w in G to be the minimum number of edges on a path from v 
to w in G; if there is no such path, we define d G (v,w) = oo. A graph G = (V,E) 
is a flow network if it has two distinguished vertices, a source s and a sink t, and 
a positive real-valued capacity u(v,w) for each edge (v,w) £ E. We extend the 
capacity function to all vertex pairs by defining u (v, w) = if (v, w) g E. A flow f 
on G is a real-valued function on vertex pairs satisfying the following constraints: 

f(v, w) < u(v, w) for all (v, w) G V x V (capacity constraint), (1.1) 

f(v,w) = -f(w,v) for all (v,w) G V x V (antisymmetry constraint), (1.2) 



^2 f(v, w) = for all v 6 V - {s, t} (flow conservation constraint). (1.3) 
wev 



Remark: The antisymmetry constraint (1.2), which is nonstandard, has two pur- 
poses: (i) it eliminates the possibility of having positive flow on both edges of an 
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opposing pair (v,w) and (w,v), a possibility that creates certain technical diffi- 
culties, and (ii) it simplifies the formal expression of constraints such as the flow 
conservation constraint (1.3). To gain an intuition one should think only of the 
positive part of the flow function; the appropriate interpretation of the flow conser- 
vation constraint is that the total flow into any vertex v £ {s,t} equals the total 
flow out of v. 

The value \f\ of a flow / is the net flow into the sink: 

l/l = £/(M). 
vev 

A maximum flow is a flow of maximum value. 

The problem we wish to solve is that of computing a maximum flow in a given 
network. Our algorithm solves this problem by manipulating a preflow f on the 
network. A preflow is a real- valued function on vertex pairs satisfying (1.1) and 
(1.2) above, as well as the following weakened form of (1.3): 

]C /(*» v ) ^ ° for all u e V — {s} (nonnegativity constraint). (1.4) 

iev 

That is, the total flow into any vertex v / s is at least as great as the total flow out 
of v. We define the flow excess e(v) of a vertex v to be J2iev /(*>*>), the net flow 
into v. 

The preflow algorithm works by examining vertices other than s and t with 
positive flow excess and pushing excess from them to vertices estimated to be closer 
to the sink t, with the goal of getting as much excess as possible to t. If the 
sink is not reachable from a vertex with a positive excess, however, the algorithm 
pushes this excess to vertices estimated to be closer to the source s. Eventually the 
algorithm reaches a state in which all vertices other than s and t have zero excess. 
At this point the preflow / is a flow; in fact, / is a maximum flow. 

Before describing the algorithm, we first address two issues: how to move flow 
excess from one vertex to another, and how to estimate the distance from a vertex 
to s or to t. 
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To deal with the first issue, we define the residual capacity rj(v,w) of a ver- 
tex pair (v,w) to be u(v,w) - f(v,w). If vertex v has positive excess and pair 
(v, w) has positive residual capacity, then an amount of flow excess up to 8 = 
min(e(v), r f (v,w)) can be moved from v to w by adding 8 to f(v,w) (and sub- 
tracting 8 from f(w,v)). Observe that there are two ways a pair (v,w) can have 
positive residual capacity: either (v,w) is an edge with flow less than its capacity 
(edge (v,w) is said to be unsaturated), or (w,v) is an edge with positive flow. In 
the former case, moving excess from v to w increases the flow on edge (v, w); in the 
latter case, it decreases the flow on (w,v). We call a pair (v,w) a residual edge if 
r f (v,w) > 0; the residual graph Gj = (V,Ej) for a preflow / is the graph whose 
vertex set is V and whose edge set Ej is the set of residual edges. 

The second issue is how to estimate the distance from a vertex to s or to t. For 
this purpose we define a valid labeling d to be a function from the vertices to the 

nonnegative integers and infinity 1 , such that d(s) = n, d{t) = 0, and d(v) < d(w) + l 

for every residual edge (v,w). The intent is that if d(v) < n, then d(v) is a lower 

bound on the actual distance from v to t in the residual graph G/, and if d(v) > n, 

then d(v) — n is a lower bound on the actual distance to s in the residual graph. It 

can be proved by induction that in the latter case t is not reachable from v in Gj.) 

To describe the algorithm, we also need the following definition. We call a vertex 
v active if v G V - {s,t}, d(v) < oo, and e(v) > 0. 

The maximum flow algorithm begins with the preflow / that is equal to the 
edge capacity on each edge leaving the source and zero on all other edges, and with 
some initial sink labeling d. The algorithm then repetitively performs, in any order, 
the basic operations, push and relabel, described in Figure 1.1. When there are no 
active vertices, the algorithm terminates. A summary of the algorithm appears in 
Figure 1.2. 



Our definition allows distance labels to be infinite. We shall show, however, that the labels stay 
finite throughout the execution of the algorithm. Infinite labels are introduced only to simplify the 
exposition. 
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Push(v, w). 

Applicability: v is active, r f (v,w) > and d(v) = d{w) + 1. 

Action: Send 6 = min(e(u), rj(v, u>)) units of flow from v to w as follows: 

/(»» «>) «- /(», w) + tf; /(«;, w) ♦- /(w, v) - 6; 

e(v) <— e(v) — 6] e(w) <— e(w) + ^. 



i?e/a6e/(u). 

Applicability: v is active and V-u; G F, r f (v,w) > => d(w) < rf(u;). 

Action: rf(w) «- min(„ iU , )eB/ (rf(u;) + 1). 

(If this minimum is over an empty set, d(v) <— oo.) 



Figure 1.1: Push and relabel operations. 

The basic operations modify the preflow / and the labeling d. A push from v to 
w increases f(v,w) and e(w) by 6 = min(e(v),rf(v,w)), and decreases f(w,v) and 
e{v) by the same amount. The push is saturating if r/(u, iw) = after the push and 
nonsaturating otherwise. A relabeling of u sets the label of v to the largest value 
allowed by the valid labeling constraints. 

Lemma 1.2.1 If f is a preflow, d is any valid labeling for f, and v is any active 
vertex, then either a push or a relabel operation is applicable to v. 

Proof: For any residual edge (v,w), the definition of a valid labeling implies that 
d(v) < d(w) + l. If a push is not applicable to v, then d(v) < d(w) + l for all residual 
edges (v, w). By the integrality of valid labelings, d(y) < d(w) for all residual edges 
(v,w), and a relabeling is applicable to v. | 

There is one part of the algorithm we have not yet specified: the choice of an 
initial labeling d. The simplest choice is d(s) = n and d{v) = for v G V - 
{s}. A more accurate choice (indeed, the most accurate possible choice) is d(v) = 
mm(d Gj (v,t), d Gj (v, s) + n) for v G V, where / is the initial preflow. The latter 
labeling can be computed in 0(m) time using backward breadth-first searches from 
the sink and from the source in the residual graph. The resource bounds we shall 
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Procedure Max-Flow (V,E,s,t,c)', 

{( initialization )) 

(( initialize preflow )) 

V(», w)e(V- {s}) x(V- {s}) do begin 

/(t>,u;)<-0; /(«;,«) «-0; 
end; 

Vu; 6 V do f(s,w) <— u(s,w)', 
{{ initialize labeling }) 
d(s) <— n; 
Vu G V - {s} do d(v) <- 0; 

(( "OOP )) 

while 3 a basic operation that applies do 
select a basic operation and apply it; 
return(/); 

end. 



Figure 1.2: The generic maximum flow algorithm. The running time of the algorithm 
depends on the order in which basic operations are applied and on details of the imple- 
mentation. 

derive for the algorithm are correct for any valid initial labeling. To simplify the 
proofs, we assume that the algorithm starts with the simple labeling. 



1.3 Correctness and Termination 

We shall prove that the generic algorithm is correct assuming that it terminates 
and then prove termination. 

Lemma 1.3.1 The algorithm maintains the invariant that d is a valid labeling. 

Proof: We use induction on the number of pushing and relabeling operations. The 
simple labeling used initially is valid because labels of all vertices other than s are 
zero, and all edges leaving s are saturated. Given that d is a valid labeling, a 
relabeling operation changing d{v) must produce a new valid labeling. Consider a 
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pushing operation that sends flow from v to w. This operation may add (w, v) to 
G f and may delete (v,w) from G f . Since d(w) = d(v) - 1, the addition of (w,v) 
to G f does not affect the invariant that d is a valid labeling. The deletion of (v, w) 
removes the corresponding constraint, which also leaves the labeling valid. | 

To prove correctness, we use the concept of an augmenting path. An augmenting 
path is a simple path from s to tin the residual graph G f . Our proof of correctness 
is based on the following classical theorem of Ford and Fulkerson [20]: 

Theorem 1.3.2 A flow f is maximum if and only if there is no augmenting path, 
i.e. t is not reachable from s in Gf. 

Lemma 1.3.3 If f is a preflow and d is any valid labeling for f , then the sink t is 
not reachable from the source s in the residual graph Gf. 

Proof: Assume by way of contradiction that there is an augmenting path s = 
u ,t;i...,t>, = t. Then / < n and (vi,v i+1 ) G E f for < i < I. Since d is a 
valid labeling, we have d(v t ) < d(v i+1 ) + 1 for < i < I. Therefore, we have 
d(s) < d(t) + 1 < n, since d(t) = 0, which contradicts d(s) = n. | 

Theorem 1.3.4 Suppose that the algorithm terminates and all distance labels are 
finite at termination. Then the preflow f is a maximum flow. 

Proof: If the algorithm terminates and all distance labels are finite, all vertices in 
V — {s,t} must have zero excess, because there are no active vertices. Therefore / 
must be a flow. This flow is maximum by Lemma 1.3.3 and Theorem 1.3.2. | 

Now we show that the algorithm terminates and that the distance labels stay 
finite during the execution of the algorithm. First we prove the following lemma: 

Lemma 1.3.5 If f is a preflow and v is a vertex with positive excess, then the 
source s is reachable from v in the residual graph Gf. 
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Proof: Let S be the set of vertices reachable from v in Gj, and suppose s g S. Let 
S = V — S. The choice of S implies that for every vertex pair i, w with i £ S and 
w 6 S, we have f{w,i) < 0. Thus 

Eves c(») = J2wev,ies / («>, *') 

= E we s,i 6 s /Hi + Eu,,igs /(«>, 

< 0. 

The term £i>es /( w 5 m tne second line equals zero by antisymmetry. Since / is 
a preflow, we have e(i) = for all i £ S, and in particular, we have e(v) = 0. | 

Lemma 1.3.6 For any vertex v, the distance label d(v) never decreases. An appli- 
cation of a relabeling operation to v increases d(v). 

Proof: Since the labeling d is changed using relabeling operations only, it is enough 
to prove the second statement of the lemma. Suppose a relabeling operation is 
applicable to v. Then for all w such that (v,w) G Ej, we have d(w) > d(v), which 
implies that min(„ iW ) € £; / (d(ti;) + 1) > d(u), so the relabeling must increase d(v). | 

We have shown that an application of the relabeling operation to a vertex in- 
creases the vertex label. The next lemma shows that the labels cannot increase 
too much. In particular, the lemma implies that the labels stay finite during an 
execution of the algorithm. 

Lemma 1.3.7 At any time during the execution of the algorithm and for any vertex 
v e V, d(v) < 2n- 1. 

Proof: The lemma is trivial for v = s and v = t. Suppose v 6 V - {s, t}. Since the 
algorithm changes only labels of active vertices, it is enough to prove the lemma for 
an active vertex v. If v is active, then e(v) > 0, so by Lemma 1.3.5 there is a simple 
path from t; to s in G f . Let v = v , ui, . . . , v t = s be such a path. The length / of 
the path is at most n - 1. Since d is a valid labeling and (u,-,u,- +1 ) e Ej, we have 
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d(vi) < d(v i+1 ) + 1. Therefore, we have d(v) = d(v ) < d(v,) + I < d(s) + (n - 1) = 
2n-l. | 

Lemma 1.3.7 allows us to amortize the work done by the algorithm over increases 
in vertex labels. The next two lemmas bound the number of relabelings and the 
number of saturating pushes. 

Lemma 1.3.8 The number of relabeling operations is at most 2n — l per vertex and 
at most (2n — l)(n - 2) < 2n 2 overall. 

Proof: Relabeling operations apply only to vertices v E V — {s,t}. A relabeling of 
v increases d(v). The label d(v) is zero initially, and the label can grow to at most 
2n - 1. Therefore there are at most 2n - 1 relabelings of each vertex in V - {s,t}, 
and the total number of relabelings in at most (2n — l)(n — 2). | 

Lemma 1.3.9 The number of saturating push operations is at most 2nm. 

Proof: For any pair of vertices v and w, consider the saturating pushes from v to w 
and from w to v. If there are any such pushes, it must be the case that (v, w) £ E 
or (w, v) G E. Consider a saturating push from v to w. In order to push flow from 
v to w again, the algorithm must first push flow from w to v, which cannot happen 
until d(w) increases by at least two. Similarly, d(v) must increase by at least two 
between saturating pushes from w to v. Since d(v) + d(w) > 1 when the first push 
between v and w occurs and d(v) + d(w) < An - 3 when the last such push occurs 
(by Lemma 1.3.7), the total number of saturating pushes between v and w is at 
most 2n — 1. Thus the total number of saturating pushes is at most 2n — 1 per edge, 
for a total over all edges of at most (2n — l)m < 2nm. | 

Next we bound the number of nonsaturating pushing operations. 
Lemma 1.3.10 The number of nonsaturating pushing operations is at most An 2 m. 
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Proof: Let <£ = Z){„|„i s active} d(v). Each nonsaturating push from a vertex v to 
another vertex w causes $ to decrease by at least one, since the push makes v 
inactive and d(w) = d(v) — 1. A saturating pushing operation causes $ to increase 
by at most 2n — 1. The total increase in $ due to saturating pushes is at most 
(2n - 1) x 2nm by Lemma 1.3.9. The total increase in $ over the entire algorithm 
due to relabeling operations is at most (2ra — l)(n-2) by Lemma 1.3.8. Immediately 
after the initialization $ is zero, $ is always nonnegative, and at the end of the 
algorithm $ is zero. Thus the total decrease in $, and hence the total number 
of nonsaturating pushing operations, is equal to the total increase in $, which is 
at most (2n - l)2ram + (2n - l)(ra - 2) < 4n 2 m for n > 4 (recall the assumption 
m > n — 1). | 

Theorem 1.3.11 The generic algorithm terminates after 0(n 2 m) basic operations. 

Proof: Immediate from Lemmas 1.3.8, 1.3.9, and 1.3.10. | 

The running time of the generic algorithm depends upon the order in which 
basic operations are applied and the details of implementation, but it is clear that 
any reasonable sequential implementation of the algorithm will run in polynomial 
time. In the next section we discuss one possible implementation with an 0(n 2 m) 
time bound. We also show that a particular ordering of the basic operations yields 
an 0(n 3 ) time bound. 

We conclude this section with a discussion of a variant of the generic maximum 
flow algorithm. Let us recall a classical concept from network flow theory, that of a 
cut. A cut S, S is a partition of the vertex set V (that is, SU ~S = V and Spl^ = 0) 
such that s G S and t G S. The capacity of the cut is 

u(S,S)= 5Z u(v,w). 

A cut is minimum if it has minimum possible capacity. The max- flow, min-cut 
theorem of Ford and Fulkerson [20,19] states that the value of a maximum flow 
equals the capacity of a minimum cut. 
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In many applications in which the maximum flow problem occurs, only the 
maximum flow value or a minimum cut is needed, not an actual maximum flow [61]. 
For such applications our maximum flow algorithm can be modified to compute a 
minimum cut and the maximum flow value without actually computing a maximum 
flow. What we have to say about the maximum flow algorithm in the remainder 
of the chapter applies to this cut algorithm as well. The only change necessary is 
to redefine an active vertex to be a vertex v G V - {s,t} such that e(v) > and 
d(v) < n. When the modified algorithm terminates, the excess e(t) at the sink is 
the value of a maximum flow, and the cut S,~S such that 5 contains exactly those 
vertices from which t is reachable in G 5 is a minimum cut [37]. For this variant of 
the algorithm, the bounds in Lemmas 1.3.8-1.3.10 can be improved by roughly a 
factor of two. 



1.4 Sequential Implementation 

Any reasonable implementation of the generic maximum flow algorithm runs in 
polynomial time. Some implementations, however, are more efficient than others. 
We shall start with a simple implementation and then refine it to improve efficiency. 

As a first step toward obtaining an efficient sequential implementation, we shall 
describe a simple refinement of the generic algorithm that runs in 0(n 2 m) time, 
matching Dinic's bound. We need some data structures to represent the network 
and the preflow. We call an unordered pair {v,w} such that (v,w) G E or (w,v) G 
E an undirected edge of G. We associate the three values u(v,w), u(w,v), and 
f(v,w)(= —f(w,v)) with each undirected edge {v,w}. Each vertex v has a list of 
the incident undirected edges {v,w}, in fixed but arbitrary order. Thus each edge 
{v,w} appears in exactly two lists, the one for v and the one for w. Each vertex v 
has a current edge {v,w}, which is the current candidate for a pushing operation 
from v. Initially, the current edge of v is the first edge on the edge list of v. The 
refined algorithm repeats the push/relabel operation, described in Figure 1.3, until 
there are no active vertices. (We shall discuss the maintenance of active vertices 
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Push/Relabel(v) . 

Applicability: v is active. 

Action: Let {v, w} be the current edge of v. 

If push(v, w) is applicable then push(v, w) 
else 

if {v, w} is not the last edge on the edge list of v then 

replace {v, w} as the current edge of v by the next edge on the edge list of v; 

else begin 

make the first edge on the edge list of v the current edge; 
relabel(v); 
end. 



Figure 1.3: The push/relabel operation. 



later). 



The push/relabel operation combines the basic operations using the above data 
structures. When applied to an active vertex v, the operation tries to push excess 
along the current edge (u, w), or advance the current edge if the pushing operation 
does not apply to the current edge. Advancing the current edge is impossible if this 
edge is the last edge on the edge list of v. In this case, pu$h/relabel makes the first 
edge on the edge list of v the current edge and apply the relabeling operation to v. 

We need to show that push/relabel uses the relabeling operation correctly. 

Lemma 1.4.1 The push/relabel operation does a relabeling only when the relabeling 
operation is applicable. 

Proof: Push/relabel applies the relabeling operation at a vertex v only when v 
is active. Just before the relabeling, for each edge (v,w), either d(v) < d(w) or 
r f (v,w) = 0, because the distance label d(v) has not changed since (v,w) was 
the current edge of v, d(w) never decreases, and rj(v,w) cannot increase unless 
d(w) > d(v). The lemma follows from the definition of a relabeling operation. | 

The refined algorithm needs one additional data structure, a set Q containing 
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Discharge. 
Applicability: Q ^ 0. 

Action: Remove the vertex v on the front of Q. 

(Vertex v must be active.) 
Repeat 

push/relabel(v); 

if w becomes active during this push/relabel operation then 
add w to the rear of Q; 
until e(v) = or d(v) increases. 
If v is still active then add v to the rear of Q. 



Figure 1.4: The discharge operation. 

all active vertices. Initially Q = {w E V - {s,t}\u(s,w) > 0}. Maintaining Q takes 
only 0(1) time per push/relabel operation. (Such an operation applied to an edge 
{v,w} may require adding w to Q and/or deleting v.) 

Theorem 1.4.2 The refined algorithm runs in 0{nm) time plus 0(1) time per 
nons Obturating pushing step, for a total of 0{n 2 m) time. 

Proof: Let v be a vertex in V - {s,t}, and let A„ be the number of edges on the 
edge list of v. Relabeling v requires a single scan of the edge list of v. By Lemma 
1.3.8, the total number of passes through the edge list of v is at most An - 1, one 
for each of the at most (2n — 1) relabelings of u, one before each relabeling as 
the current edge runs through the list, and one after the last relabeling. Every 
push/relabel operation selecting v either causes a push, changes the current edge 
of v, or increases d{v). The total time spent in push/relabel operations selecting v 
is 0(nA„) plus 0(1) per push out of v. Summing over all vertices and applying 
Lemmas 1.3.9 and 1.3.10 gives the theorem. | 

To obtain a better running time we need to reduce the number of nonsaturating 
pushes. We do this in a way similar to that used by Shiloach and Vishkin [66]. 
Namely, we exploit the freedom we have in selecting vertices for push/relabel oper- 
ations by using a first-in, first-out selection strategy; i.e., we maintain Qasa queue. 
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The first-in, first-out algorithm consists of applying the discharge operation until 

Q is empty. The discharge operation consists of applying push/relabel operations 

to an active vertex at least until the excess becomes zero or the label of the vertex 
increases. 

There is still some flexibility in this algorithm, namely in how long we keep 
applying push/relabel operations to a vertex v. Figure 1.4 describes one extreme 
case, in which we stop as soon as e(v) = or v is relabeled. At the other extreme we 
can continue until v becomes inactive, which may involve several relabelings of v. 
In the sequential case, our analysis is valid for both extremes and all intermediate 
variants. In the parallel case, our analysis applies only to the version described in 
Figure 1.4, but it can be easily modified to handle other cases as well. 

To analyze the first-in, first-out algorithm, we need to introduce the concept 
of passes over the queue. Pass one consists of the discharging operations applied 
to the vertices added to the queue during the initialization. Given that pass i is 
defined, pass i + 1 consists of the discharging operations applied to vertices on the 
queue that were added during pass i. 

Lemma 1.4.3 The number of passes over the queue is at most 4n 2 . 

Proof: Define the potential function $ to be $ = ma,x{d(v)\v is active}. Consider 
the effect on $ of a single pass over the queue. If no distance label changes during 
the pass, each vertex succeeds in moving its excess to lower-labeled vertices, so $ 
decreases during the pass. If $ is not changed by the pass, some vertex label must 
increase by at least one. If $ increases, some vertex label must increase by at least 
as much as $ increases. The total number of passes in which $ stays the same or 
increases is thus at most 2n 2 by Lemma 1.3.7. Since $ = initially and $ is always 
nonnegative, the total number of passes in which $ decreases is also at most 2n 2 . 
Hence the total number of passes is at most 4n 2 . | 



Corollary 1.4.4 The number of nonsaturating pushes during the first-in, first-out 
algorithm is at most 4n 3 . 



1.5. USE OF DYNAMIC TREES 27 

Proof: There is at most one nonsaturating push per vertex in V — {s,t} per pass. 
I 

Theorem 1.4.5 The first-in, first-out algorithm runs in 0(n 3 ) time. 
Proof: Immediate from Theorem 1.4.2 and Corollary 1.4.4. | 

An alternative strategy for vertex selection, which we call the maximum distance 
method, is to always select a vertex v in Q with d(v) maximum. This strategy also 
gives an 0(n 3 ) running time, as a proof similar to that of Theorem 1.4.5 shows. 

1.5 Use of Dynamic Trees 

We have now matched the 0(n 3 ) time bound of Karzanov's algorithm. To obtain a 
better bound, we must reduce the time per nonsaturating pushing operation below 
0(1). We do this by using the dynamic tree data structure of Sleator and Tarjan 
[68,69,71]. This data structure allows us to maintain a set of vertex-disjoint rooted 
trees in which each vertex v has an associated real value g(v ), possibly oo or — oo. 
We regard a tree edge as directed toward the root, i.e. from child to parent. We 
denote the parent of a vertex v by p(v). We adopt the convention that every vertex 
is both an ancestor and a descendant of itself. The tree operations we shall need 
are described in Figure 1.5. 

The total time for a sequence of / tree operations starting with a collection of 
single- vertex trees is 0(7 log k), where k is an upper bound on the maximum number 
of vertices in a tree. (The implementation of dynamic trees presented in [69,71] does 
not support find-size operations, but it is easily modified to do so. See [37].) 

In our application the edges of the dynamic trees are a subset of the current 
edges of the vertices. The current edge {v, w} of a vertex v e V - {s,t} is eligible to 
be a dynamic tree edge (with p(v) = w) if d(v) = d{w) + 1 and r f (v, w) > 0. Not all 
eligible edges are tree edges, however. The value g(v) of a vertex v in its dynamic 
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find-root(v): Find and return the root of the tree containing vertex v. 

find-size(v): Find and return the number of vertices in the tree containing 
vertex v. 

find-value(v): Compute and return g(v). 

find-min(v): Find and return the ancestor w of v of minimum value g(w). 
In case of a tie, choose the vertex w closest to the root. 

change-value(v,x): Add real number x to g{w) for each ancestor w of v. (We 
adopt the convention that oo + (— oo) = 0.) 

link(v,w): Combine the trees containing vertices v and w by making w 
the parent of v. This operation does nothing if v and w are 
in the same tree or if v is not a tree root. 

cut(v): Break the tree containing v into two trees by deleting the edge 

from v to its parent. This operation does nothing if v is a tree 
root. 

Figure 1.5: Dynamic tree operations. 

tree is rj(v,p(v)) if v has a parent and oo if v is a tree root. Initially, each vertex 
is in a one- vertex dynamic tree and has value oo. We limit the maximum tree size 
to k, where k is a parameter to be chosen later. 

By using appropriate tree operations we can push flow along an entire path in 
a tree, either causing a saturating push or moving flow excess from some vertex 
in the tree all the way to the tree root. By combining this idea with a careful 
analysis, we are able to show that the number of times a vertex is added to Q is 
0(nm + n 3 /k). At a cost of O(log k) for each tree operation, the total running time 
of the algorithm is 0((nm + n s /k)logk), which is minimized to within a constant 
factor at 0(nmlog(n 2 /m)) for the choice k — n 2 /m. 

The details of the improved algorithm, which we call the dynamic tree algorithm, 
are as follows. The heart of the algorithm is the procedure send(v) defined in Figure 
1.6, which pushes excess from a nonroot vertex v to the root of its tree, cuts edges 
saturated by the push, and repeats these steps until e(v) = or v is a tree root. 

At the top level, the dynamic tree algorithm is exactly the same as the first-in, 
first-out algorithm of Section 1.4: we maintain a queue Q of active vertices and 
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Send(v). 

Applicability: v is active. 

Action: While find-root(v) ^ v and e(v) > do begin 

send 6 <- mm(e(v), find-value(find-min(v))) units of flow 
along the tree path from v by performing change-value(v , —8); 
while find-value(find-min(v)) = do begin 
i <— find-min(v); 

perform cut(i) followed by change-value(i, oo); 
end; 
end. 



Figure 1.6: The Send operation. 

repeatedly perform discharging operations until Q is empty. However, we replace 
the push/relabel operation with the tree-push/relabel operation described in Figure 
1.7. 

A tree -push/relabel operation applies to an active vertex v that is the root of 
a dynamic tree. There are two main cases. The first case occurs if the current 
edge {v, w} of v is eligible for a pushing operation. If the trees containing v and w 
together have at most k vertices, we link these trees by making w the parent of v 
and do a send operation from v. If these trees together contain more then k vertices, 
we do an ordinary pushing operation from v to w followed by a send from w. The 
second case occurs if the edge {v, w} is not eligible for a pushing operation. In this 
case we update the current edge of v and relabel v if necessary. If v is relabeled, we 
cut all tree edges entering u, to maintain the invariant that all dynamic tree edges 
are eligible for pushing operations. 



It is important to realize that this algorithm stores values of the preflow / in two 
different ways. If {v, w} is an edge that is not a dynamic tree edge, f(v, w) is stored 
explicitly, with {v,w}. If {v,w} is a dynamic tree edge, with w the parent of v, 
then g(v) = u(v, w) - f(v, w) is stored implicitly in the dynamic tree data structure. 
Whenever a tree edge (v, w) is cut, g(v) must be computed and f(v, w) updated to 
its current value. In addition, when the algorithm terminates, flow values must be 
computed for all edges remaining in dynamic trees. 
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Tree-Push/Relabel(v) . 

Applicability: v is an active tree root. 

Action: Let {v, w} be the current edge of v. 

(1) If d{v) = d(w) + 1 and r f (v, w) > then begin 
(la) If find-size(v) + find-size(w) < k then begin 

make w the parent of v by performing 

change-value(v, -oo), change-value(v , r f (v , w)), and link(v,w)\ 
push excess from v to w by performing send(v); 
end 
(lb) else (( find-size(v) + find-size(w) > k)) begin 

apply a pushing operation to move excess from v to w; 
perform send(w); 
end 
end; 

(2) else {(d(v) < d(w) or r f (v,w) = 0)) 

(2a) if {v, w} is not the last edge on the edge list of v then 

replace {v, w} as the current edge by the next edge on the list 
(2b) else (({v,w} is the last edge on the edge list of v)) begin 

make the first edge on the list the current one; 

perform cut(i) and change-value(i) for every child i of v; 

apply a relabeling operation to v; 
end. 



Figure 1.7: The tree-push/relabel operation. 

Two observations imply that the dynamic tree algorithm is correct. First, any 
edge {v,w} that is in a dynamic tree has d(v) = d(w) + 1. Therefore in case (la) 
of tree-push/relabel, vertices v and w are in different trees, and the algorithm never 
attempts to link a dynamic tree to itself. Second, a vertex v that is not a tree 
root can have positive excess only in the middle of case (1) of a tree-push/relabel 
operation. To see this, note that only in this case does the algorithm add excess to 
a nonroot vertex, and this addition of excess is followed by a send operation that 
moves the nonroot excess to one or more roots. 

Lemma 1.5.1 The dynamic tree algorithm runs in 0(nmlogk) time plus O(logib) 
time per addition of a vertex to Q. 

Proof: The condition in subcase (la) of tree-push/relabel guarantees that the max- 
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imum size of any dynamic tree is k. Thus the time per dynamic tree operation 
is 0(logA;). Each tree-push/relabel operation takes 0(1) time plus 0(1) tree oper- 
ations plus 0(1) tree operations per cut operation (in invocations of send and in 
subcase (2b)) plus time for relabeling (in subcase (2b)). The total relabeling time 
is 0{nm). The total number of cut operations is at most the number of link opera- 
tions, which is at most 2nm by a proof like that of Lemma 1.3.9. (Another way of 
getting a bound on the number of cut and link operations is to observe that a cut 
operation corresponds to a saturating push or to an edge scan during relabeling, 
and the number of link operations exceeds the number of cut operations by at most 
n - 1.) The total number of tree-push/relabel operations is 0(nm) plus one per 
addition of a vertex to Q. Combining these observations gives the lemma. | 

We define passes over the queue Q exactly as in Section 1.4. The proof of Lemma 
1.4.3 remains valid, which means that the number of passes is at most An 2 . 

The next lemma is the crucial part of the analysis. 

Lemma 1.5.2 The number of times a vertex is added to Q is 0(nm + n 3 /k). 

Proof: A vertex v can be added to Q only after d(v) increases, which happens 
at most 2n times, or as a result of e{v) increasing from zero, which can happen 
only in subcases (la) and (lb) of tree-push/relabel. In either subcase, the number 
of vertices added to Q is at most one more than the number of cuts performed 
during the invocation of send in the subcase. Thus the number of additions to 
Q in subcases (la) and (lb) is at most 2nm (the maximum number of cuts) plus 
the number of occurrences of the subcases. There are at most 2nm occurences of 
(la) (the maximum number of links). There are at most 2nm occurrences of (lb) 
in which the invocation of send(w) causes a cut, and at most 2nm occurrences of 
(lb) in which the push from v to w is saturating. Let us call an occurrence of (lb) 
nonsaturating if it adds a vertex to Q but causes neither a cut nor a saturating 
push. It remains for us to count the number of nonsaturating occurrences. 

We need a few definitions. For any vertex u, we denote the dynamic tree con- 
taining u by T u and the number of vertices it contains by \T U \. Tree T u is small if 
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\T U \ < k/2 and large otherwise. At any time, and in particular at the beginning of 
any pass, there are at most 2n/k large trees. 

Consider a nonsaturating occurrence of (lb) during a given pass, say pass i. The 
condition in (lb) guarantees that either T v or T w is large, giving us two cases to 
consider. 

Suppose T v is large. Vertex v is the root of T v . The nonsaturating occurrence 
of (lb) removes all the excess from u, which means that a nonsaturating occurrence 
can apply to a given vertex v only once during a given pass. If T v has changed since 
the beginning of pass i, we charge the occurrence of (lb) to the link or cut that 
changed T v most recently before the occurrence. The number of such occurrences 
over all passes is at most one per link and two per cut, for a total of at most 6nm. 
(A link forms one new tree; a cut, two.) If T v has not changed since the beginning 
of pass i, we charge the occurrence of (lb) to T v . Since T v is large and there are 
at most 2n/k large trees at the beginning of pass i, there are at most 2n/k such 
charges per pass, for a total of at most 4n 3 /k over all passes. 

Suppose on the other hand that T w is large. The occurrence of (lb) adds the 
root of T w , say r, to Q (otherwise this occurrence of (lb) need not be counted). 
A given vertex r can be added to Q at most once during a given pass. If T w has 
changed since the beginning of pass i, we charge the occurrence of (lb) to the link 
or cut that changed T w most recently before the occurrence. The number of such 
occurrences over all passes is at most Qnm. If T w has not changed, we charge the 
occurrence to T w . The number of such charges over all passes is at most 4n 3 /k. 

Summing our estimates, we find that there are at most 2n 2 + 20nm + 8n 3 /k 
additions to Q altogether, giving the lemma. | 

Theorem 1.5.3 The dynamic tree algorithm runs in 0(nm\og(n 2 j 'm)) time if k 
is chosen equal to n 2 /m. 

Proof: Immediate from Lemmas 1.5.1 and 1.5.2. | 
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As in Section 1.4, we can replace first-in, first-out selection of vertices for dis- 
charging steps by maximum distance selection, and still obtain the same running 
time bound. 

1.6 Parallel and Distributed Implementation 

The synchronous parallel version of our algorithm is a modification of the first-in, 
first-out algorithm of Section 1.4. The algorithm proceeds in pulses, each of which 
consists of a number of operations applied in parallel. Each pulse is divided into 
three stages. The pushing of flow is done during the first stage, the relabeling of 
vertices is done in the second stage, and flow pushed to a vertex in the first stage 
is added to its excess in the third stage. We make three changes in the algorithm. 
First, we restrict the algorithm so that it stops processing a vertex v as soon as 
e(v) = or v is relabeled. Second, instead of using a queue for selection of vertices 
to be processed, we process all active vertices in parallel. Third, the flow pushed to 
a vertex v during a parallel step is not added to e(v) until the third stage. To be 
more precise, the parallel version consists of repeating the ■pulse operation described 
in Figure 1.8 until there are no active vertices. 

The parallel algorithm is almost a special case of the first-in, first-out algorithm, 
the only difference being in the values used in relabeling and flow excess compu- 
tations: in the first-in, first-out algorithm, these computations in pass i use the 
most recent label and excess values, some of which may have been computed earlier 
in pass i. Nevertheless, a proof just like that of Lemma 1.4.3 gives the following 
analogous result for the parallel algorithm: 

Lemma 1.6.1 The number of pulses made by the parallel algorithm is at most An 2 . 

Corollary 1.6.2 The number of nonsaturating pushes made by the parallel algo- 
rithm is at most 4n 3 . 

For the distributed implementation of this algorithm, our computing model is 
as follows [31,2]. We allow each vertex v of the graph to have a processor with an 
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Procedure Pulse. 

For all active vertices v in parallel do begin 
(( stage 1)) 

push flow from v until e(v) = or \/w such that d(w) — d(v) — 1, r/(v,w) = 0. 
« stage 2)) 
If e(v) > then begin 

d'(v) «- mm w]r} ( ViW)>0 (d(w) + 1); 
if d(v) ^ d'(v) then begin 
d(v) «- d'(v)i 

broadcast d(v) to all neighbors of v; 
end; 
end. 

({ stage 3)) 

Add flow pushed to v in stage 1 to e(v). 
end. 



Figure 1.8: The Pulse operation. 

amount of memory proportional to A„, the number of neighbors of v. This processor 
can communicate directly with the processors at all neighboring vertices. We assume 
that local computation is much faster than inter-processor communication. Thus 
as a measure of computation time we use the number of rounds of message-passing. 
We are also interested in the total number of messages sent. 

A synchronous distributed implementation of the parallel algorithm works as 
follows. Each vertex processed during a pulse sends updated flow values to the 
appropriate neighbors. New vertex labels are also transmitted to neighbors, but 
only at the end of the pulse. Since flow always travels in the direction from larger 
to smaller labels, this delaying of the label broadcasting guarantees that flow only 
travels through an edge in one direction during a pulse. An easy analysis shows that 
in the synchronous case the distributed algorithm takes 0(n 2 ) rounds of message- 
passing and a total of 0(n 3 ) messages. 

For parallel implementation, our computing model is a PRAM [21] without 
concurrent writing. The implementation in this model is very similar to that of 
the distributed implementation, except that computations on binary trees must be 
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performed to allow each vertex to access its incident edges fast. Because of these 
binary trees, each pulse takes 0(log n) time, and the parallel time of the algorithm 
is O(n 2 \ogn). The ideas of Shiloach and Vishkin [66] apply to our algorithm to 
show that 0(n) processors suffice to obtain the 0(n 2 log n) time bound. See [66] for 
details. We describe a different parallel implementation of the algorithm in Chapter 
3. 

Now we discuss two implementations of the algorithm in the asynchronous dis- 
tributed model of parallel computation [31,2]. Awerbuch (private communication, 
1985) has observed that in the asynchronous case, the synchronization protocol of 
[2] can be used to implement our algorithm in 0{n 2 log n) rounds and 0(n 3 ) mes- 
sages. The same bounds can be obtained for the Shiloach- Vishkin algorithm [66], 
but only by allowing more memory per processor: the processor at a vertex v needs 
0(nA v ) storage. Vishkin (private communication) has reduced the space required 
by this algorithm to a total of 0(n 2 ) (from 0{nm) ). Nevertheless, our algorithm 
has an advantage in situations where memory is at a premium. 

Our algorithm can be modified to work in the asynchronous model without the 
use of the synchronization protocol, achieving better running time but using more 
messages. This asynchronous version of the algorithm synchronizes locally using 
acknowledgments. When a vertex v pushes its excess to vertex w such that, accord- 
ing to the local information at v, d(v) = d(w) + 1, it sends a message (v, 8, d(v)) and 
updates e(i>). The vertex v will not push flow to w again or change d(v) until v re- 
ceives an acknowledgment from w. When a vertex w receives a message (u, 8, d(v)), 
it first checks if d{y) = d(w) + 1 (because the value of d(w) in the processor v at 
the time it sent the message may be out of date). If d(v) — d(w) + 1, then w sends 
to v a message of the form (accept, w, 6,d(w)). Otherwise it sends to v a message 
of the form (reject, w, 8, d(w)), where d(w) is the correct value of the distance label 
of w. The accepting or rejecting messages serve as acknowledgments. In addition, 
a rejecting message causes v to update its excess, its local value of d(w), and d(v) 
if necessary. When a distance label of a vertex increases, it informs its neighbors 
about the new value of the label. 
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Theorem 1.6.3 The asynchronous distributed implementation of the algorithm that 

uses acknowledgments runs in 0(n 2 ) time using 0(n 2 m) messages and 0(A„) mem- 
ory per processor. 

Proof: To analyze the message complexity of the algorithm, note that the total 
number of messages is the number of messages generated by the distance label in- 
creases plus twice the number of (accepting or rejecting) acknowledge messages. 
The number of messages generated by the distance label increases is at most 2nm. 
There is at most one rejecting message per edge per distance label increase, for a to- 
tal of 2nm (by Lemma 1.3.7). The same arguments as in the proofs of Lemmas 1.3.9 
and 1.3.10 give Inm and 4ro 2 m bounds on the number of accepting messages cor- 
responding to saturating and nonsaturating pushes. The total message complexity 
of the algorithm is thus 0(n 2 m). 

To bound the running time of the algorithm, we need to introduce a unit of 
time. Given an execution of an algorithm, we define a time unit to be the longest 
time from the point when a message is originated by a sender to the point when the 
message is processed by the receiver. For example, a push- acknowledgment pair of 
operations takes two units of time. Note that if during a time interval (t, t + 4] no 
vertex label increases, then the $ function defined as in the proof of Lemma 1.4.3 
must decrease during this time interval. To see this, observe that during the time 
interval [t + l,i + 4) each vertex has the correct information about the distance 
labels of its neighbors, so all pushes initiated during the time interval [t + l,t + 2) 
are accepted by time t + 4. A proof similar to that of Lemma 1.4.3 yields an 0(n 2 ) 
bound on the running time of the algorithm. | 

An alternative way to obtain a fast distributed or parallel algorithm is to 
use a parallel version of maximum distance selection: during each pulse, apply 
push/relabel steps to every active vertex v for which d{y) is maximum. This re- 
quires a preprocessing step at the beginning of each pulse to compute the maximum 
d(v), but it simplifies other calculations, since during a given pulse a vertex cannot 
both send and receive flow, which allows the computations of flow excess to proceed 
concurrently with the push/relabel steps. 
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1.7 Remarks 



Our concluding remarks concern three issues: (i) better bounds, (ii) exact distance 
labeling, and (iii) efficient practical implementation. Regarding the possibility of 
obtaining better bounds for the maximum flow problem, it is interesting to note 
that the bottleneck in the sequential version of our algorithm is the nonsaturat- 
ing pushes, whereas the bottleneck in the parallel version is the saturating pushes. 
Recently, Ahuja and Orlin [1] have devised a scaling algorithm based on the ap- 
proach described in this chapter. Their algorithm runs in 0(nm + n 2 log U) time 
(assuming that the edge capacities are integers not exceeding U). This improves 
Gabow's bound of 0(nm log U) mentioned in the introduction. We wonder whether 
an 0(nm) sequential time bound can be obtained through more careful handling of 
the nonsaturating pushes, possibly avoiding the use of the dynamic tree data struc- 
ture. Perhaps also an 0(m(logn) k ) parallel time bound can be obtained through 
the use of a parallel version of the dynamic tree data structure. 

It is possible to modify our algorithm so that when a pushing operation is 
executed, each distance label is exactly the distance to the sink or to the source 
in the residual graph (i.e. if d(v) < n, then d(v) — do f (v,t); if d(v) > n then 
d(v) = dG f (v,s) + n). The modification involves a stronger interpretation of the 
current edge of a vertex, which should be unsaturated and lead to a vertex with 
a smaller label. If a pushing step saturates the current edge of v, a new current 
edge is found by scanning the edge list of v and relabeling if the end of the list is 
reached, as in a push/relabel step. If a relabeling step changes the label of v, the 
current edge must be updated for each vertex i such that (i,v) is the current edge 
of i. One can show that these computations take 0(nm) time in total during the 
algorithm. 

Whether maintaining exact distance labels improves the practical performance 
of the algorithm is not clear, because the work of maintaining the exact labels may 
exceed the extra work due to nonexact labels. The above observation, however, 
suggests that as long as we are interested in an fl(nm) upper bound on an imple- 
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mentation of the generic algorithm, we can assume that the exact labeling is given 
to us for free. 

Our algorithm is practical. In chapter 3, we describe an implementation of 
the algorithm and experimental results obtained using the implementation. The 
algorithm was also used in the context of computational physics [58]. In a prac- 
tical implementation it is important to make the algorithm as fast as possible. 
We offer a heuristic that may speed up the algorithm. The heuristic periodically 
updates the distance labels by performing breadth-first searches backwards from 
the sink and source in the residual graph. These searches compute, for each ver- 
tex v, the distances do f (v,s) and d,G f (v,t). The new distance label of v is set to 
min(c?G / (u, s), do f (v, t) + n). There are several possible strategies for deciding when 
to recompute labels. One is to do so after every n relabeling operations. Another is 
to do so every time an edge into the sink is saturated or an edge out of the source 
has its flow reduced to zero. Neither of these strategies affects the asymptotic time 
bound of the algorithm, but they may improve its practical performance. 

Another important issue in a practical implementation is what strategy to use 
for selecting vertices for discharging steps. Although the best theoretical bounds 
we have obtained are for first-in, first-out and maximum distance selection, other 
strategies, such as last-in, first-out and maximum excess, deserve consideration. 
Ahuja-Orlin scaling algorithm [1], for example, selects vertices with excess larger 
than a certain threshold h in such a way that no excess greater than 2h is created 
during a scaling iteration. 



Chapter 2 



The Minimum- Cost Flow Problem 



2.1 Introduction 

In this chapter we introduce a method that allows to apply techniques developed 
for the maximum flow problem to the more general minimum-cost flow problem. 
The method incorporates the idea of cost scaling used in the algorithms of Rock [62] 
and of Bland and Jensen [8], and the concept of relaxed complementary slackness 
conditions embodied in the algorithms of Bertsekas [7,6] and of Tardos [70]. We 
show how to use both the maximum flow approach of Dinic [14] and the maximum 
flow approach of Chapter 1 in the context of this method. 

The minimum-cost flow problem is that of finding a feasible flow of minimum 
cost in a network with capacity constraints and costs on edges. This problem has 
a wider range of applications then the maximum flow problem discussed in the 
previous chapter. Extensive discussion of the problem and its applications appear 
in the books of Ford and Fulkerson [19], Lawler [52], Papadimitriou and Steiglitz 
[60], and Tarjan [71]. 

All known polynomial-time algorithms for the problem are based on the idea of 
scaling. 1 This idea was introduced by Edmonds and Karp [15], who used the idea 
to design the first polynomial-time algorithm for the problem. Scaling algorithms 



1 For applications of scaling to other problems, see [26]. 
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for the minimum-cost flow problem work by solving a sequence Po,Pi, . . . ,Pl of 
minimum-cost flow problems. The problem P, is obtained by considering the i most 
significant bits of the capacities (capacity scaling) or the costs (cost scaling). The 
first problem Pq is easy to solve because all the capacities or costs are equal to zero. 
A solution of the problem P, helps to solve the next problem P,-+i. If all of the 
relevant values have at most L bits, the solution of Pl is the desired solution. 

Table 2.1 summarizes polynomial-time algorithms for the minimum-cost flow 
problem. The running time of the algorithms is given in terms of the number n of 
vertices, the number m of edges, the maximum absolute value U of capacities, and 
the maximum absolute value C of costs C. When U (or C) appears in a bound, 
the capacities (or the costs) are assumed to be integer. When stating the running 
times, we assume the best time bounds known for the shortest path and maximum 
flow subroutines used by some of these algorithms: 0(m + nlogn) for the shortest 
path subroutine [22] and 0(nmlog(n 2 /m)) for the maximum flow subroutine (see 
Chapter 1). Algorithms 1, 2, 5, 6, and 8 use the capacity scaling and algorithms 3, 
4, 7, and 9 use cost scaling. The algorithms 4, 5, 6, and 8 are strongly polynomial: 
their running time does not depend on U or C. 

Since the running times of the algorithms in Table 2.1 are expressed in terms 
of different parameters, the algorithms cannot be compared directly. The previous 
algorithms 1-8, in the three comparable groups, rank as follows. Algorithms 1 
and 2 give the best capacity-dependent bound, algorithms 3 and 7 give the best 
cost-dependent bound, and algorithm 8 gives the best strongly polynomial bound. 
For most applications, however, one can assume that U = n°^ and C = n°^. 
Under these assumptions, the bound of O(mlog(n)(m + nlogn)) achieved by the 
capacity-scaling algorithms 1 and 2 is the best among the previous algorithms. 

In this chapter we present a general approach to the minimum-cost flow prob- 
lem. The approach combines methods for solving the maximum flow problem with 
successive approximation techniques. We use this approach to construct algorithms 
for the problem with upper bounds of 0(nmlog(n)log(nC)), 0(n 5 / 3 m 2 / 3 log(nC)), 
and 0(n 3 log(nC)), which significantly improves the best previous cost-dependent 
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# 


Date 


Discoverer 


Running Time 


References 


1 


1972 


Edmonds and Karp 


0(mlog(U)(m + nlogn)) 


[15] 


2 


1980 


Rock 


0(m\og(U)(m + nlogn)) 


[62] 


3 


1980 


Rock 


0(n log(C)(nm log(n 2 /m))) 


[62] 


4 


1984 


Tardos 


0(m 4 ) 


[70] 


5 


1984 


Orlin 


0(m 2 log(n)(m + nlogn)) 


[59] 


6 


1985 


Fujishige 


O(m 2 log(n)(m + nlogn)) 


[23] 


7 


1985 


Bland and Jensen 


0(nlog(C)(nmlog(n 2 /m))) 


[8] 


8 


1986 


Galil and Tardos 


O(n 2 log(n)(m + nlogn)) 


[30] 


9 


1987 


Goldberg and Tarjan 


0(min(nm log n, n 5 / 3 m 2 / 3 , n 3 ) log( nC)) 


[38] 



Table 2.1: Polynomial-time algorithms for the minimum-cost flow problem. Algorith 
9 is presented in this chapter. 



m 



bound achieved by algorithms 3 and 7 in the table, as well as the best bound under 
the assumptions U = n°^ and C = n°^ discussed above. The algorithm gives the 
best bound on the complexity of the problem for C — o(n n ) and C = o(U n /n). Our 
approach is in many respects similar to the approaches taken by Bland and Jensen 
[8] to develop algorithm 7, and by Bertsekas [6] to develop an exponential-time al- 
gorithm for the problem. Our techniques are more powerful, however, and lead to 
more efficient algorithms. 

The new algorithm works by successive approximation. It starts by finding an 
approximate solution and then iteratively improves the current solution, each time 
doubling the quality of approximation by halving the error parameter e. The inner 
loop subroutine that improves the approximation is based on generalizations of 
techniques for solving the maximum flow problem. When the error parameter is 
small enough, the current solution is optimal, and the algorithm terminates. To 
measure the quality of a solution, we use the notion of e-optimality, which is related 
to the classical technique of perturbing a linear programming problem to avoid 
degeneracy (see, for example, [32]). The notion of e-optimality is motivated by 
the relaxation of the complementary slackness conditions [7,70]. The termination 
condition used in our algorithm is due to Bertsekas [6]. 
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Our approach can be viewed as a generalization of the cost-scaling approaches 
of Rock [62] and of Bland and Jensen [8]. The optimal solution to the problem P, 
solved by a cost-scaling algorithm is an e-optimal solution to the original problem 
with e = 2' logC l +1_ *, so the error decreases by a factor of two from one scaling 
iteration to another. In fact, earlier versions of our algorithm used the cost-scaling 
approach. However, the use of true costs throughout the algorithm simplifies its 
analysis and implementation. 

Our successive approximation approach is different from the traditional scaling 
approaches in one important way. In traditional scaling, each problem P, is solved 
exactly, even though the data used to define the problem P, is imprecise. Our 
approach can be interpreted as solving the imprecise subproblems approximately, 
with the error parameter being within a constant factor of the precision of the 
problem data. The work saved by solving the intermediate problems approximately 
is one of the reasons for the improved efficiency of our method. This observation 
can be used to improve other scaling algorithms. The minimum-cost flow problem is 
closely related to other flow- like problems and to the linear programming problem, 
so the techniques developed in this chapter may apply to the other problems as 
well. 

Some of the results of this section are based on observations of Robert Tarjan. 
These results include Lemma 2.6.5, which simplifies implementations of the generic 
subroutine and the corresponding analysis, and a way to combine the successive 
approximation approach with Dinic's layered network approach, which leads to the 
results described in section 2.9. The results presented in this chapter will also 
appear in [38]. 



2.2 Definitions and Notation 

In this section we define the minimum- cost circulation problem and introduce the 
notation and terminology used throughout the chapter. The minimum-cost circula- 
tion problem is a generalization of the maximum flow problem discussed in Chapter 
1. The minimum-cost flow problem is also a special case of the linear program- 
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ming problem and is usually defined in linear programming terms. Although we 
use several theorems which have their roots in the theory of linear programming, 
most arguments presented in this chapter are graph theoretic. Consequently, we 
formulate the problem in graph-theoretic terms. Our formulation of the problem 
is equivalent to other formulations of the minimum-cost flow and minimum-cost 
circulation problems that can be found in the books and papers referred in the 
introduction to this chapter. 

A circulation network is a directed graph G = (V, E) with upper and lower 
capacity bounds and costs on edges. Circulation networks are different from flow 
networks defined in Section 1.2 in three ways: circulation networks have costs on 
edges, lower bounds on edge capacities, but they do not have sources or sinks. As 
before, we denote the size of V by n and the size of E by m, and we assume that 
m > n — 1 > 4 (as in Section 1.2). We call an unordered pair {v,w} such that 
(v,w) G E or (w,v) G E an undirected edge of G. For notational convenience, we 
extend the capacity functions and the cost function to all pairs of vertices. Let R 
denote the set of real numbers. The capacity bounds are given by functions u : 
V x V — ► R and I : V x V — > R with the following constraints for all (v, w) eV xV: 

l(v,w) < u(v,w) (consistency constraints), (2-1) 

u(v,w) = —l(w,v) (capacity antisymmetry constraints), (2-2) 

l(v, w) = u(u, w) = if (u, w) $ E and (w, v) # E. (2.3) 

Property (2.3) is required to extend the capacity functions to all pairs of vertices. 

A pseudoflow 2 is a function f : V x V ^> R satisfying the following constraints 
for all (u, w) G V x V: 

l(»,to) < f(v,w) < u(v,w), (capacity constraints), (2-4) 



2 The concept of a pseudoflow is different from the preflow concept of Karzanov defined in Section 
1.1 in that the flow conservation constraints are completely dropped. 
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f(v,w) = —f(w,v) (flow antisymmetry constraints). (2-5) 

A circulation is a pseudoflow that satisfies 



^/(u,u;) = (conservation constraints) (2-6) 



for all v G V . 

We assume that the costs of edges are given by a cost function c : V x V — > R 
that satisfies the following constraints for all (v,w) GVxV: 

c(v,w) = —c(w,v) (cost antisymmetry constraints). (2-7) 

We extend the cost function to pairs of vertices by defining c(v, w) = for all (v , w) 
such that (v,w) £" E and (w,v) £ E. The cost of a circulation f is given by the 
following expression: 



\ £ c(v,w)f(v,w). (2.8) 

(v,w)ev*v 



(The factor of 1/2 appears because we count the cost of the flow between each pair 
of vertices twice.) The minimum- cost circulation problem is to find a circulation of 
minimum cost (an optimal circulation). 

Remark: We refer to equations (2.2), (2.5), and (2.7) as the antisymmetry con- 
straints. One should think of a positive and a negative direction for each undirected 
edge of G, with the capacity and cost constraints given for the positive direction 
and derived for the negative direction using the antisymmetry constraints. 

Another important concept related to the problem is the concept of vertex prices. 
To gain intuition, consider a vertex v and a real number x. Suppose that we add 
x to prices of all edges going into v and subtract x from prices of all edges going 
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out of v. Because of the conservation constraints at v, the cost of / is not changed 
by this transformation (the cost of each unit of flow going into v increases by x, 
and the cost of each unit of flow going out of v decreases by x). Therefore the 
transformed problem is equivalent to the original one. To define the prices formally, 
we introduce a price function p : V — * R, and define the price of a vertex v to be 
p(v). The reduced cost function c p is defined by c p (v,w) = c(v,w) — p(v) + p(w). 
In the linear programming interpretation of the problem, the prices correspond to 
dual variables. 

Given a pseudoflow /, we define the residual capacity function rf : V x V — » R 
by rf(v,w) — u(v,w) — f(v,w). The residual graph Gf = (V,Ef) is the directed 
graph with vertex set V containing all edges with positive residual capacity: Ef = 
{(v,w)\rf(v,w) > 0}. The balance bf(v) of a vertex v, is the difference between the 
incoming and outgoing flows, or, more formally, the function 6/ : V x V — »■ R defined 
by bf(v) = "%2 W £V f( w i v )- If / is a circulation, then bf(v) — for all v. Given a 
pseudoflow /, we say that a vertex v is active if bj{y) > 0. Note that J2vev K v ) = 
for any pseudoflow, so a pseudoflow is a circulation if and only if there are no active 
vertices. 

We also need the following standard definitions. An augmenting path (cycle) is 
a simple path (cycle) in Gf. The cost of a path (cycle) is the sum of the costs of all 
edges on the path (cycle). 

2.3 Optimality and Approximate Optimality 

In this section we define the notion of an e-optimal pseudoflow and show that for 
e < 1/ro, an e-optimal circulation is optimal. 

The following theorem of Ford and Fulkerson [19] provides an optimality crite- 
rion for a circulation. 

Theorem 2.3.1 A circulation f is a minimum-cost circulation if and only if there 
exists a price function p such that V(v,w) £VxV, 
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c p (v, w) > =► f(v, w) = l(v, w) (2.9) 



ant 



c p (v,w) <0=> f(v,w) = u(v,w). (2.10) 

The optimality conditions (2.9) and (2.10) are called complementary slackness 
conditions, and an edge (u, w) satisfying these conditions is said to be in kilter. 
We use the term kilter because of the relationship between our algorithm and the 
out-of-kilter method [57,24]. Figure 2.1a shows a kilter diagram [52], which is a 
pictorial representation of the complementary slackness conditions. 

The antisymmetry constraints (2.2), (2.5), and (2.7) imply that an edge (v,w) 
is in kilter if and only if the corresponding edge (w, v) is in kilter. 

The following theorem [19] gives an optimality criterion that does not involve 
the price function. 

Theorem 2.3.2 A circulation f is optimal if and only if the residual graph Gf 
contains no cycles of negative cost. 

To define e-optimality, we use the notion of e-relaxation of the complementary 
slackness conditions [7,70]. Given e > 0, we say that a pseudoflow / is e-optimal if 
the following relaxations of the complementary slackness conditions hold: V(u, w) E 
VxV, 

c p (v,w) > e =>• f(v,w) = l(v,w) (2-11) 



and 



c p (v,w) < — e =$> f(v,w) = u(v,w). (2.12) 
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(a) Kilter diagram 



(b) Approximate optimality 



Figure 2.1: 

(a) Kilter diagram. If / is an optimal circulation, the edge (u, w) is on the heavy curve. 

(b) e-optimality. If / is e-optimal, the edge (v,w) can be in the shaded rectangle as 
well as on the heavy curve. 

We say that a pseudoflow / is e-optimal if there exists a price function p with 
respect to which the pseudoflow / is e-optimal. Figure 2.1b illustrates the concept 
of e-optimality in terms of kilter diagrams. 

The antisymmetry constraints (2.2), (2.5), and (2.7) imply the following lemma: 

Lemma 2.3.3 An edge (v,w) satisfies conditions (2.11) and (2.12) if and only if 
the corresponding edge (w,v) satisfies conditions (2.11) and (2.12). 



The following simple fact is used extensively in our presentation. 



48 CHAPTER 2. THE MINIMUM-COST FLOW PROBLEM 

Lemma 2.3.4 Suppose that pseudoflow f is e-optimal with respect to a price func- 
tion p. Then for any residual edge (v,w) we have c p (v,w) > — e (i.e., the cost of a 
residual edge cannot be too small). 

The following theorem of Bertsekas [6] shows that when e is small enough, an 
e-optimal circulation is optimal. 

Theorem 2.3.5 Assume that all edge costs are integers. Then for any < e < 1/n, 
an e-optimal circulation f is optimal. 

Proof: Consider a cycle in Gf. By Lemma 2.3.4, the e-optimality of / implies that 
the cost of the cycle is at least — ne, which is greater than — 1. Since the costs are 
integers, the cost of the cycle must be at least 0. The theorem then follows from 
the optimality criterion given by Theorem 2.3.2. | 



2.4 High-Level Description of the Algorithm 

Theorem 2.3.5 suggests the algorithm Min-Cost summarized in Figure 2.2. First, 
the algorithm finds a circulation (using a maximum flow algorithm) and sets the 
prices of all vertices to zero. The resulting circulation is C-optimal (recall that C is 
the largest absolute value of an edge cost). Then, the algorithm iteratively improves 
the approximation (using the Improve- Approximation subroutine), until the error 
becomes less than 1/n. At this point, the current solution is optimal. 

Remark: The algorithm need not use a maximum flow subroutine in the initializa- 
tion stage. It can start with any pseudoflow. If the problem is infeasible, we can 
discover this fact during the first execution of Improve- Approximation because the 
increase in vertex prices will be greater than allowed by the analysis below. We use 
the maximum flow subroutine only to be able to assume, without loss of generality, 
that the input problem is feasible, so that we do not have to worry about feasibility 
in our presentation. 
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Procedure Min-Cost(V,E,l,u,c); 

e *- max( ViW ) eE \c(v,w)\i 

f e <— feasible circulation; (( use a maximum flow subroutine )) 

Vi>, Pt(v) <- 0; 

while e > 1/n do 

(A/25 Pc/2) <— Improve-Approximation(f t ,p e ,e)', 

e «- € /2; 
end; 
return(/ e ); 

end. 



Figure 2.2: The minimum-cost flow algorithm. 

Theorem 2.4.1 Let D(n,m) be the running time of the Improve- Approximation 
subroutine. Then the minimum-cost flow algorithm runs in 0(D(n,m)\og(nC)) 
time and returns a minimum-cost flow. 

Proof: Immediate from Theorem 2.3.5 and the above discussion. | 

The inputs to the Improve- Approximation subroutine are a circulation f € , a price 
function p £ , and an error parameter e, such that f t is e-optimal with respect to p t . 
The outputs of the subroutine are a circulation f t / 2 and a price function p t / 2 such 
that / £ / 2 is (e/2)-optimal with respect to p £ / 2 . The subroutine sets the flow through 
every out-of- kilter edge to the upper or lower bound to bring the edge in kilter. 
The resulting pseudoflow is (e/2)-optimal (in fact, 0-optimal), but the conservation 
constraints (2.6) may be violated. Then, the pseudoflow is transformed into a 
circulation by applying a sequence of operations that preserves (e/2)-optimality. 
This transformation is done using generalizations of maximum flow techniques. As 
we shall see, both the techniques described in Chapter 1 and the techniques of Dinic 
[14] can be used. 

The cost scaling algorithms of Rock and of Bland and Jensen use a maximum 
flow algorithm in the inner loop. Since the Improve- Approximation subroutine is a 
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generalization of a maximum flow algorithm, our minimum-cost flow algorithm is 
similar to the cost scaling algorithms. The cost-scaling algorithms halve the error 
in 0(n) iterations of the inner loop, however, whereas our algorithm halves the the 
error in approximation in a single iteration. 

Intuitively, the algorithm can be viewed, as simulated annealing [48]. At the 
beginning of each iteration, we release the potential energy of the system, and then 
let the system cool with the speed determined by e. When e is big, the energy is 
dissipated quickly, but we are unlikely to find an optimal solution. When e is small 
enough, we are guaranteed to find the optimal solution. At each iteration, the value 
of e is selected so that the system cools reasonably quickly. 



2.5 Generic Improve- Approximation Subroutine 

The generic Improve- Approximation subroutine is a generalization of the generic 
maximum-flow algorithm of Section 1.2. An implementation of the generic subrou- 
tine using the discharge operation as described in Section 2.7 is almost identical 
to the earlier minimum-cost flow algorithm of Bertsekas [6], but we use the sub- 
routine in a different way. In our algorithm the subroutine is used to reduce the 
approximation parameter e by a factor of two at each iteration. In the algorithm 
of Bertsekas, the value of e is set to l/(n + 1) at the initialization, so the algorithm 
terminates in one iteration — but may take exponential time. 

Our presentation and analysis of the generic subroutine and its implementations 
follow the presentation and analysis of the generic maximum flow algorithm. Al- 
though the generalization of the algorithms is straight-forward, the generalization 
of the analysis is more complicated. In fact, some time bounds that we prove in 
this chapter are not quite as good as the corresponding bounds for the maximum 
flow algorithms. The differences are discussed in Section 2.7. 

The generic Improve- Approximation subroutine is described in Figure 2.3. The 
subroutine forces all edges in kilter and then transforms the resulting pseudoflow 
into an (e/2)-optimal circulation. The basic operations used to manipulate the 
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Procedure Improve- Approximation^ f,p, e); 

(( initialization)) 

V(t>, w) £ V X V do begin 

if c p (v, w) > then /(u, w) <— /(u, w); 

if c p (v, w) < then /(u, u;) <— m(u, u>); 
end; 

(( "°0P)> 

while 3 a basic operation that applies 

select a basic operation and apply it; 

return(/,p)j 



end. 



Figure 2.3: The generic Improve-Approximation subroutine. The running time of the 
subroutine depends on the order in which basic operations are selected and on details 
of implementation. 

pseudoflow and the prices are push and update-price. The operations are de- 
scribed in Figure 2.4 and are illustrated in terms of kilter diagrams in Figure 
2.5. A push through an edge (v,w) increases f(v,w) and bj(w) by at most 8 = 
min(bf(v),rf(v, w)) and decreases f(w,v) and bj(v) by the same amount. For 
the purpose of the analysis, we distinguish between saturating and nonsaturat- 
ing pushes. A push is saturating if rj(v,w) = after the push and nonsaturating 
otherwise. The update-price operation sets the price of a vertex v to the highest 
value allowed by the (e/2)-optimality constraints. 

The following lemmas give important properties of the push and update-price 
operations. 

Lemma 2.5.1 A pushing operation preserves the (e/2)-optimality of a pseudoflow. 

Proof: Let / be an (e/2)-optimal pseudoflow with respect to a price function p. Sup- 
pose a pushing operation is applied to an edge (v, w). Since (v, w) £ Ej and a push- 
ing operation is applicable to (v,w), it follows that — e/2 < c p {v,w) < — e/4. Note 



52 CHAPTER 2. THE MINIMUM-COST FLOW PROBLEM 



Push{y, w). 

Applicability: v is active, rf(v,w) > 0, and c p (v,w) < —e/4. 

Action: Send 6 — min(6/(v), rj(v, w)) units of flow from v to w. 



Update- Price (v) . 

Applicability: v is active and Vty £ V rj(v,w) > =>■ c p (v,w) > —e/4. 

Action: Replace p(v) by mm( v ^ eE (p(w) + c(v, w) + e/2). 



Figure 2.4: Push and update-price operations. 

that the relaxed complementary slackness conditions do not restrict flow through 
edges of small cost, therefore changing the flow through any residual edge (v,w) 
with \c p (v, w)\ < e/2 preserves (e/2)-optimality. It follows that a pushing operation 
preserves (e/2)-optimality. | 

Lemma 2.5.2 Suppose f is an (e/2) -optimal pseudoflow with respect to a price 
function p and a price updating operation is applied to a vertex v. Then the price 
of v increases by at least e/4 and the pseudoflow f is (e/2) -optimal with respect to 
the resulting price function p'. 

Proof: First we prove that the price of v increases by at least e/4. The price 
updating operation changes the price of the vertex v from p(v) to p'(v) and does 
not affect prices of other vertices. Since a price updating operation is applicable to 
u, we have c(v, w) — p(v) + p(w) > —e/4 for all residual edges (v, w). Therefore, we 
have c(v,w) + p(w) + e/2 > p(v) + e/4 for all edges (v,w) G Ef, so the new price 
p'(v) > p(v) + e/4 by definition of the price updating operation. 

Now we prove the second claim of the lemma. Because of the antisymmetry, 
it is enough to show that for all w £ V, the relaxed complementary slackness 
conditions corresponding to (v,w) remain valid. The constraint 2.11 holds after 
the price update because it holds before the price update and from the first part 
of the proof we know that c p i(v,w) < c p (v,w). The constraint 2.12 clearly holds 
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if f(v,w) = u(v,w). If f(v,w) < u(v,w) then we have (v,w) E Ef so c(v,w) — 
p'(v) + p'(w) = c(t>, w) — p'(v) + p(w) > — e/2 by the definition of the price updating 
operation. Therefore the constraint 2.12 also holds in this case. | 

Lemma 2.5.3 If a pseudoflow f is (e/2)-optimal with respect to a price function 
p and v is an active vertex, then either a pushing or a price updating operation is 
applicable to v. 

Proof: Suppose v is active and no pushing operation is applicable to v. By (e/2)- 
optimality of /, for any residual edge (v,w), we have c p (v,w) > —e/2 by Lemma 
2.3.4. Furthermore, since the edge (v, w) cannot be used for pushing, we have 
c p (v,w) > — e/4. Therefore a price updating operation is applicable. | 

The generic version of Improve- Approximation, described in Figure 2.3, repeti- 
tively performs an applicable basic operation on the current pseudoflow / and price 
function p until no basic operation applies. We shall prove that the generic ver- 
sion of Improve- Approximation is correct if it terminates, and then we shall prove 
termination. 

Theorem 2.5.4 // the generic Improve- Approximation subroutine terminates on 
an e-optimal input circulation, then the pseudoflow f output by the subroutine is an 
(e/2)-optimal circulation. 

Proof: The subroutine terminates when there are no vertices with positive balance 
and therefore no vertices with negative balance, so / is a circulation. The (e/2)- 
optimality of / follows from Lemmas 2.5.1 and 2.5.2 and the fact that the pseudoflow 
computed during the initialization step of Improve- Approximation is (e/2)-optimal 
(in fact, 0-optimal). | 

We conclude this section with a remark that shows how to find an optimal 
circulation given an optimal price function (a price function p is optimal if there is 
a circulation / that satisfies the complementary slackness conditions 2.9 and 2.10). 
Suppose that we force the edges in kilter in the same way as in the initialization 
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c p i v > w )h 



s/z 



-e/2 



1(v,vjI 



y 




u(y,w} 



f(v,v)f 



Figure 2.5: Push(i;,u>) is allowed if b f (v) > and (v,w) is in the shaded rectangle 
(including the heavy curve). If b f (v) > but no push operation is applicable, Update- 
Price^) is applicable. Increasing the price of v corresponds to shifting the kilter diagram 
down. The absence of edges in the shaded rectangle guarantees that the price of v 
increases by at least e/4. 

step of the generic Improve- Approximation subroutine. It can be shown using the 
techniques of the next section that the resulting pseudoflow can be converted into 
a circulation by changing the flow only through edges with zero reduced cost. The 
resulting circulation satisfies the complementary slackness conditions and therefore 
is optimal. This circulation can be found in one maximum flow computation. 



2.6 Analysis of the Generic Subroutine 



In this section we analyze the generic Improve- Approximation subroutine. The 
analysis is similar to the analysis of the generic maximum flow algorithm. We 
start the running time analysis by bounding the amount of increase in vertex prices 
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during an execution of the generic subroutine. 

Lemma 2.6.1 Let f be a pseudoflow and let f be a circulation. Then for any 
vertex v with a positive balance bj(v), there exists a vertex w with a negative balance 
bj(w) and a sequence of vertices u\, . . . , u;_i such that (u, wj, . . . , ttj_i, w) is a simple 
path in Gj and (u;, u;_i, . . . , u\, v) is a simple path in Gf>. 

Proof: Fix a vertex v with a positive balance. Let G + — (V,E+) where E+ = 
{(ij)lf'(ij) > f(hj)} and let G_ = (V,E.) where E_ = {(i,j)\f'(i,j) < f(ij)}. 
This definition implies that E + C Ej, because for any edge (i,j) G Ef we have 
f(i,j) < f'(i,j) < u(i,j). Similarly, we have 25_ C Ej>. Furthermore, if (i,j) is an 
edge in G+, then (j,i) is an edge in G- by antisymmetry. Therefore it is enough to 
show that if 6/(u) > 0, then there is a simple path (u,«i, . . . , ui-i, w) in G + such 
that bf(w) < 0. 

Assume by the way of contradiction that such a path does not exist, i.e. all 
vertices reachable from v in G + have nonnegative balance. Let S be the set of 
all vertices reachable from v in G+ and let S = V — S. The vertices in S have 
nonnegative balance, and since v e S and b/(v) > 0, we have Yljes KJ) > 0- Also, 
for every pair of vertices (i,j) G S x S, we have f(i,j) > f'(i,j), for otherwise we 
would have j G S. 

Now we obtain a contradiction as follows. 

= E(,j)esxs/'(*>i) (/' is a circulation) 

< E ( ,j )e sxs/(*.i) 

= £(,,j) 6 sxs/(*>j) + £(«\j)esxs/(*,j) (antisymmetry) 

= E(tj)esxv/(*,j) 

= —JZjesbf(j) (definition of balance) 

< 

I 

The following key lemma bounds the amount of increase in prices during an 
execution of the subroutine. This lemma is similar to the lemma of Bland and 
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Jensen [8] that bounds the number of maximum flow computations in a scaling 
step. 

Lemma 2.6.2 The price of any vertex v increases by at most (3/2)ne during an 
execution of Improve- Approximation. 

Proof: Since we change vertex prices only by using price updating operations, which 
apply only to vertices with positive balance, it is enough to prove the lemma for a 
vertex v that has a positive balance at some time during the execution of Improve- 
Approximation. Let / and p be the (e/2)-optimal pseudoflow and price function at 
this point. Recall that f e and p t are the e-optimal circulation and the price function 
at the beginning of the execution of Improve- Approximation. Let v, u\, . . . , tx/_i, w 
be a sequence of vertices satisfying Lemma 2.6.1. Applying the e-optimality condi- 
tions to the edges on the path Pj € = (w, u;_i, . . . , iti, v) in Gf c , we obtain 



p t (w)<p e (v) + h + J2 <hJ)- ( 2 - 13 ) 

(tj) on p u 



Applying the (e/2)-optimality conditions to the edges on the path 
Pf = (u,iii, . . . , itj_i,u>) in Gf, we obtain 



p(v) < pH + /(6/2) + E( j ,,-)onp / c(j,0 ( . 

= p(t0) + Z(6/2)-E W ) on p /€ c(i,j). V- l *> 



By Lemma 2.6.1, the balance b/(w) is negative. Note that application of basic 
operations cannot create a new vertex with a negative balance. Thus w has had a 
negative balance after the initialization step of Improve- Approximation. The prices 
of vertices with negative balance do not change, so p(w) = p £ (w). Combining 
this observation with inequalities 2.13 and 2.14, we get p(v) < p t (v) + (3/2)/e < 
p € (v) + (3/2)ne. I 

.Remaric: Although Lemmas 2.6.1 and 2.6.2 are sufficient for our presentation, the 
following statement gives additional insight into the problem. This statement can 
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be proven in a way similar to Lemmas 2.6.1 and 2.6.2. Let fi be an ei- optimal 
circulation with respect to a price function pi and let f% oe an 62-optimal circulation 
with respect to a price function p2. Let S be a set of vertices such that for all s £ S 
we have pi(s) — P2(s) and any vertex v E V is reachable from some vertex s £ S in 
G$ x or in Gf 2 . Then for any vertex v, we have \p\{v) — P2(v)\ < (ei + e 2 )n. 

Lemma 2.6.3 The number of the price updates during an execution of Improve- 
Approximation is at most 6n 2 . 

Proof: Immediate from lemmas 2.5.2 and 2.6.2. | 

Lemmas 2.6.2 and 2.6.3 enable us to amortize the operations performed by the 
algorithm over increases in vertex prices. 

Lemma 2.6.4 The number of the saturating push operations during an execution 
of Improve- Approximation is at most Inm. 

Proof: For any undirected edge {v, w}, consider saturating pushes from v to w and 
from w to v. Consider a saturating push from v to w. In order to push flow from 
v to w again, the subroutine must first push from w to v, which can not happen 
until the price of w increases by at least e/2. Similarly, p{y) must increase by at 
least e/2 between saturating pushes from w to v. By charging all saturating pushes 
from v to w (except for the first one) to price increases of v and applying Lemma 
2.6.2, we can bound the number of pushes that use the undirected edge {v ,w} by 
2 + 6n < 7n (assuming n > 2). Summing over all undirected edges gives the desired 
bound. | 

At this point we introduce the concepts of an admissible edge and the admissible 
graph. Given a pseudoflow / and a price function p, we say that an edge (v,w) is 
admissible if (v,w) E Ej and c p (v,w) < — e/4. Note that if (v,w) is an admissible 
edge and v is active, then a pushing operation is applicable to (v,w). For a given 
price function p, the admissible graph Ga is the graph of all admissible edges: 
G A = (V,E A ), where E A = {(v,w) E E f \c p (v,w) < -e/4}. 



58 CHAPTER 2. THE MINIMUM-COST FLOW PROBLEM 

The following lemma is due to Robert Tarjan. 

Lemma 2.6.5 At any time during the execution of the generic Improve- Approxi- 
mation subroutine, the admissible graph Ga is acyclic. 

Proof: We prove the lemma by induction on the number of basic operations. For 
the basis, note that after the initialization there are no admissible edges, because all 
edges in the residual graph have nonnegative costs. Now suppose that Ga is acyclic 
before a basic operation is applied. A pushing operation can delete an admissible 
edge but cannot create a new admissible edge. Therefore a pushing operation cannot 
create a cycle in Ga- 

Now consider a price updating operation applied to a vertex v. This operation 
affects only admissible edges adjacent to v, so it is enough to show that after the 
operation there is no cycle in Ga that passes through v. Before the price update, 
the reduced cost of any residual edge entering v is at least — e/2. The price update 
increases this cost by at least e/4. Therefore after the price update, there are no 
admissible edges going into v. It follows that after the price updating operation, 
there is no cycle in Ga that passes through v. | 

The following lemma bounds the number of nonsaturating pushes. This lemma 
has been strengthened to its current form with the help of Ron Rivest. 

Lemma 2.6.6 The number of the nonsaturating pushing operations in an execution 
of Improve- Approximation is 0{n 2 m). 

Proof: For the purpose of the proof, we define a function h that maps vertices into 
real numbers. For a vertex v, let h(v) be if the out-degree of v in the admissible 
graph is zero. Otherwise, let h(v) be the minimum, over all paths P in the admissible 
graph that start at v, of the sum of reduced costs of edges in P. Define $ = if 
there are no active vertices and $ = Yli v \ v is active) ^( u ) otherwise. Note that for 
any vertex v, we have —(e/2)n < h(v) < and therefore — (e/2)n 2 < $ < 0. Since 
the admissible graph is acyclic, each nonsaturating pushing operation causes $ to 
increase by at least e/4. 
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Next we bound the total amount of decrease in $. A saturating pushing op- 
eration causes $ to decrease by at most (e/2)n. The total decrease in $ due to 
saturating pushes is at most 7n 2 ra(e/2) by Lemma 2.6.4. As we have seen in the 
proof of Lemma 2.6.5, a price updating operation applied to a vertex v removes 
all admissible edges going into v. Therefore, such an operation can decrease h(v), 
but cannot increase h(w) for any w =fi v. It follows that each price update causes 
$ to decrease by at most (e/2)n. The total amount of decrease in $ due to price 
updates is at most 6rc 3 (e/2) by Lemma 2.6.3. After the initialization step, $ is at 
most ra 2 (e/2), and $ is always nonnegative. Therefore the number of nonsaturating 
pushes is at most 14ra 2 m + 12n 3 + 2n 2 = 0{n 2 m). | 

The following theorem bounds the number of basic operations in the generic 
Improve- Approximation subroutine. 

Theorem 2.6.7 The generic Improve- Approximation subroutine terminates after 
0(n 2 rri) basic operations. 

Proof: Immediate from Lemmas 2.5.1, 2.5.2, 2.6.3, 2.6.4, and 2.6.6. | 

As we know, the maximum flow problem is a special case of the minimum-cost 
circulation problem. The standard reduction transforms an instance of the maxi- 
mum flow problem into an instance of the minimum-cost flow problem as follows. 
Lower bounds and costs of all edges are set to zero. Then, a new edge (t,s) going 
from the sink t to the source s is introduced. The upper bound on the flow through 
this edge is set to the sum of capacities of all edges going out of the source, the 
lower bound on the flow is set to zero, and the cost of the edge is set to —n. Since 
this new edge is the only edge with negative cost and all other edges have zero cost, 
an optimal circulation maximizes flow through this edge — and maximizes flow 
from the source to the sink in the original network. If we saturate the edge (t,s), 
the resulting pseudoflow is 1-optimal with respect to the zero price function. The 
generic Improve- Approximation subroutine applied at this point works almost in the 
same way as the generic maximum flow algorithm. The only difference is that the 
subroutine changes the price of the source, whereas the maximum flow algorithm 
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does not change the distance label of the source. This difference is minor, because 
the maximum flow algorithm can be easily modified to allow changes in the distance 
label of the source [37]. 

We would like to conclude this section with a discussion of the choice of the 
value —e/4 as an upper bound on the reduced cost of admissible edges. There are 
other possible choices for this lower bound. In particular, we can set this bound 
to zero, allowing pushes through residual edges with a negative reduced cost and 
modifying the price updating operation to apply to a vertex v only when reduced 
costs of all residual edges going out of v are nonnegative. The asymptotic complexity 
bounds stated in this chapter on the performance of the generic subroutine and its 
derivatives still apply if we make this change in the algorithm (although some proofs 
and constant factors change). 

The justification for the choice of the value —e/4 is on purely intuitive level. 
We can define the cost of a pseudoflow by equation (2.8). The choice of the value 
—e/4 (or — ke for a constant k such that < k < 1/2) assures that moving a unit of 
flow by a pushing operation changes the cost of a pseudoflow by at least —e/4 (or 
— ke times the number of flow units pushed). This change in the cost of pseudoflow 
seems important for the validity of Conjecture 1 stated in the next section and for 
potential generalization of the Ahuja-Orlin algorithm [1] to the minimum-cost flow 
problem. 



2.7 Sequential Implementation 

The running time of the generic subroutine depends upon the order in which the 
basic operations are applied and on the details of implementation, but it is clear 
that any reasonable sequential implementation will run in polynomial time. As 
a first step toward obtaining an efficient sequential implementation of the generic 
subroutine, we shall describe a simple refinement of the subroutine that runs in 
0(n 2 m) time. Then, we describe a first-in, first-out version of the generic subrou- 
tine, which is similar to the first-in, first-out maximum flow algorithm described in 
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Push/update(v) . 
Applicability: v is active. 

Action: Let {v, w} be the current edge of v. 

If Push(i;, w) is applicable, apply it; 
else 

if {v, w}is not the last edge on the edge list of v then 
replace {v, w} as the current edge of v 
by the next edge on the edge list of v; 
else begin 

make the first edge on the edge list of v the current edge; 
apply Update- Price(v); 
end. 



Figure 2.6: The push/update operation. 

section 1.4. Unlike in the maximum flow case, however, we are not able to prove the 
0(n 3 ) time bound on the first-in, first-out algorithm (although we believe that this 
bound holds). We shall show a different implementation of the generic Improve- 
Approximation subroutine that runs in 0(n 3 ) time. In the next section, we shall 
describe implementation of the generic subroutine that uses the dynamic tree data 
structure and runs in 0(nm log n) time. 

We need some data structures to represent the network and the pseudoflow. We 
associate a positive direction and the following four values with each undirected edge 
{v,w}: l(v,w), u(v,w), c(v,w), and f(v,w). Each vertex v has a list of the incident 
undirected edges {v,w} in a fixed order. Each edge {v,w} appears in exactly two 
lists, the one for v and the one for w. Each vertex v has a balance bf(v) and a 
current edge {v,w}, which is the current candidate for a push out of v. Initially 
the current edge is the first edge on the edge list of v. After the initialization, the 
refined subroutine applies the push/update operation described in Figure 2.6 until 
there are no active vertices. 

We need to show that push/update uses the price update operation correctly. 
Lemma 2.7.1 The push/update procedure uses a price updating operation only 



62 CHAPTER 2. THE MINIMUM-COST FLOW PROBLEM 

when this operation is applicable. 

Proof: The push/update procedure applies the price updating operation only to 
active vertices with positive balance. Just before the price update, for each edge 
(v,w) either c p (v,w) > —e/4 or rf(v,w) — 0, because p(v) has not changed since 
{v,w} was a current edge, rj(v,w) cannot increase unless c p (v,w) > —e/4, and 
p(w) never decreases. The lemma follows from the definition of the price updating 
operation. | 

The refined subroutine needs one additional data structure, a set Q containing 
all vertices with positive balance. Initially Q contains vertices whose balance has 
been made positive during the initialization. Maintaining Q takes only 0(1) time 
per push/update operation. (Such an operation applied to an edge {u, w} may 
require adding w to Q and/or deleting v.) 

Theorem 2.7.2 The refined algorithm runs in 0(nm) time plus 0(1) time per 
nonsaturating pushing step, for a total of 0(n 2 m) time. 

Proof: Let v G V and let A„ be the number of edges on the edge list of v. A 
price update of v requires a single scan of the edge list of v. By the proof like that 
of Lemma 2.6.3, the total number of passes through the edge list of v is at most 
12n + 1: one for each of at most Qn price updates of u, one before each price update 
as the current edge runs through the list, and one after the last price update. Every 
push/update operation selecting v either causes a push, changes the current edge of 
i>, or increases p(v). The total time spent in push/update operations selecting v is 
0(nA v ) plus 0(1) time per push out of v. Summing over all vertices and applying 
Lemmas 2.6.4 and 2.6.6 gives the theorem. | 

As in the maximum flow case, our next step is to describe the first-in, first-out 
version of the generic subroutine. This version of the subroutine maintains Q as 
a queue and applies discharge operations until the queue is empty. The discharge 
operation, described in Figure 2.7, is the same as the maximum flow algorithm 
except that the push/update operation is used instead of the push/relabel operation. 
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Discharge. 
Applicability: Q ^ 0. 

Action: Remove the vertex v on the front of Q. 

(Vertex v must be active.) 
Repeat 

apply push/update(v); 

if w becomes active during this push/update operation then 
add w to the rear of Q ; 
until bf(v) = or p(v) increases. 
If v is still active then add v to the rear of Q. 



Figure 2.7: The discharge operation. 

We define passes over the queue in the same way as in Section 1.4. Pass one 
consists of the discharging operations applied to the vertices added to the queue 
during the initialization. Given that pass i is defined, pass i + 1 consists of the 
discharging operations applied to vertices on the queue that were added during 
pass i. 

Lemma 2.7.3 The number of passes over the queue in the execution of the first-in, 
first- out Improve- Approximation subroutine is 0(n 3 ). 

Proof: Define the function h in the same way as in the proof of Lemma 2.6.6, and 
define $ = if / is a circulation and $ = min^^x)} h(v) otherwise. Recall that 
for any v, we have -(e/2)n < h(v) < and therefore -(e/2)n < $ < 0. Consider 
the effect on $ of a single pass over the queue. If the price function has not changed 
during the pass, $ must increase by at least e/4. If the price function changes 
during the pass, $ decreases by at most (e/2)n. The total number of passes during 
which $ decreases is at most Qn 2 by Lemma 2.6.3, and the total amount of decrease 
in $ is at most 3n 3 e. After the initialization step $ is at least — (e/2)n, and at the 
end of the subroutine $ is zero. The number of passes during which $ increases is 
at most 12n 3 + 2n, and the overall number of passes is 6n 2 + 12n 3 + In — 0(n 3 ). 
I 
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Note that the bound given by the above lemma is not as good as the corre- 
sponding bound of 0(n 2 ) on the number of passes for the maximum flow algorithm. 
We suspect that the analysis given in Lemma 2.7.3 can be improved to match the 
corresponding maximum flow bound to within a constant factor. 

Conjecture 1 The bounds given by Lemma 2.7. S can be improved to 0(n 2 ), either 
for the described version of the first-in, first out algorithm or for a modified version. 



All implementations of the maximum flow algorithm generalize to the minimum 
cost flow case. The bounds that we are able to prove for the implementations of 
the generic Improve- Approximation subroutine that are based on the first-in, first- 
out strategy, however, are worse then the corresponding maximum flow bounds. 
Conjecture 1 would imply the same asymptotic bounds for all variations of the 
Improve- Approximation subroutine based on the first-in, first-out strategy as for 
the corresponding variations of the maximum flow algorithm. In particular, the 
conjecture would imply an 0(nralog(n 2 /m)) sequential running time bound and an 
0(n 2 ) bound on the number of parallel pulses. 

The 0(n 3 ) bound on the number of passes of the first-in, first-out algorithm is 
not strong enough to prove an 0(n 3 ) time bound on the Improve- Approximation 
subroutine. Next we show a different implementation of the subroutine that runs 
in 0(n 3 ) time. This implementation is due to Charles Leiserson. We call this 
implementation the wave subroutine because of its similarity to the maximum flow 
algorithm described in [72] . When applied to the generic maximum flow algorithm 
of Section 1.2, the wave method yields an 0(n 3 ) time algorithm as well. 

Unlike the previous implementations we have discussed, the wave implementa- 
tion does not maintain a queue of active vertices. The implementation maintains 
a list L of all vertices instead, and preserves the invariant that the vertices on the 
list are topologically ordered with respect to the admissible graph Ga (i-e., for any 
two distinct vertices v and w, if w is reachable from v in Ga, then w appears after 
v on the list). 
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Discharge- 2{v). 
Applicability: v is active. 
Action: Repeat 

Apply push/update(v); 

until bf(v) = or p(v) increases. 

If p(v) has increased then move v to the beginning of L; 



Figure 2.8: The discharge-2 operation. 

Initially, the list L contains the vertices of G in arbitrary order. The implemen- 
tation makes passes over the list, applying the discharge-2 operation, described in 
Figure 2.8, to active vertices. The discharge-2 operation consists of applying the 
push/relabel operations to an active vertex until the excess becomes zero or the price 
of the vertex increases. In the latter case, the vertex is moved to the beginning of 
the list L, but the processing of the list L continues from the previous position of 
this vertex. The subroutine terminates when there are no active vertices. 

The key to the analysis of the wave subroutine is the observation that, because 
of the topological ordering of vertices on the list L, if no price update occurs during 
a pass over the vertex list, the subroutine terminates. 

Lemma 2.7.4 The number of passes over the vertex list in the execution of the 
wave implementation of the Improve- Approximation subroutine is 0(n 2 ). 

Proof: First we show, by induction on the number of basic operations, that the 
wave subroutine maintains a topological ordering of the vertices with respect to the 
admissible graph G A (except in the middle of the discharge-2 operation). The basis 
is trivial, because immediately after the initialization stage of the generic subrou- 
tine there are no admissible edges. A pushing operation preserves the topological 
ordering, because this operation cannot create a new admissible edge. Immediately 
after the price of a vertex is changed, the vertex is moved to the beginning of the list 
L. The resulting ordering is topological because a price updating operation applied 
to a vertex v deletes all edges of Ga going into v. 
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Next we show that if the vertex prices do not change during a pass over the 
vertex list, the subroutine terminates. This implies that there is at most one pass 
during which the price function does not change. The total number of passes is at 
most 6n 2 + 1, because by Lemma 2.6.3, the number of passes during which the price 
of some vertex changes is at most 6n 2 . 

Suppose that the price function does not change during a pass. Then no price 
updating operations are performed during this pass, and therefore the ordering of 
the vertices on the list L does not change and every vertex is able to get rid of all 
of its excess. Since the vertices are processed in topological order, no vertex has 
excess after the pass, and the algorithm terminates. | 

The 0(n 2 ) bound on the number of passes allows us to prove the 0(n 3 ) bound 
on the running time of the wave subroutine. 

Theorem 2.7.5 The wave subroutine runs in 0(n 3 ) time. 

Proof: The work done by update-price and saturating push operations can be 
bounded by 0(nm) by a proof like that of Theorem 2.7.2. The number of non- 
saturating push operations is 0{n 3 ) because there are 0(n 2 ) passes and at most 
one nonsaturating push per vertex per pass. Finally, operations on the vertex list 
L require 0{n) work per pass, for a total of 0(n 3 ) work. | 



2.8 Use of Dynamic Trees 

We improve the 0(n 3 ) bound on the Improve- Approximation subroutine given by 
Theorem 2.7.5 to 0{nm log n) by using the dynamic tree data structure discussed 
in Section 1.5. The resulting subroutine is similar to the dynamic tree version of 
the maximum flow algorithm. Unlike the maximum flow algorithm, however, the 
subroutine implements a generic version that does not use the first-in, first-out 
ordering of operations, so we are not able to use the balancing of tree sizes. 
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Send(v). 

Applicability, v is active. 

Action: While find-root(v) ^ v and bj(v) > do begin 

send S = min(bf(v),find-value{find-Tnin(v))) units of flow 

along the tree path from v by performing change-value(v , — 6); 
while find-value(find-min(v)) = do begin 
u <— find-min(v); 

perform cut(u) followed by change-value(u, oo); 
end; 
end. 



Figure 2.9: The send operation. 

As in the maximum flow implementation that uses the data structure, the edges 
of the dynamic trees are a subset of the current edges of the vertices. The current 
edge {v, w} of a vertex v is eligible to be a dynamic tree edge if the edge (i>, w) is 
admissible, but not all eligible edges are tree edges. The value g(v) of a vertex v in 
its dynamic tree is rj(v,parent(v)) if v has a parent, oo if v is a tree root. Initially 
each vertex is in a one- vertex dynamic tree and has value oo. 

By using appropriate tree operations we can push flow along an entire path in 
a tree, either causing a saturating push or moving flow excess from some vertex in 
the tree all the way to the tree root. The dynamic tree operations are charged to 
the price updates and saturating pushes in such a way that each price update and 
each saturating push is charged a constant number of times. Therefore the total 
number of charges is 0{nm) and, since each dynamic tree operation costs 0(log n), 
the running time of the dynamic tree subroutine is 0(nm log n). 

The details of the improved subroutine, which we call the dynamic tree subrou- 
tine, are as follows. The heart of the subroutine is the procedure send(v) defined 
in Figure 2.9, which pushes excess from a nonroot vertex v to the root of its tree, 
cuts edges saturated by the push, and repeats these steps until bj{y) = or v is a 
tree root. 

At the top level, the dynamic tree subroutine is exactly the same as the simple 
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Tree-Push /update(y) . 

Applicability: v is an active tree root. 

Action: Let {v, w} be the current edge of v. 

(1) If edge (v,w) is admissible then begin 
make w the parent of v by performing 

change-value(v,—oo), change-value(v,rf(v,w)), and link(v,w); 
push excess from v by performing send(v); 
end; 

(2) else (( edge (v,w) is not admissible)) 

(2a) if {v, w} is not the last edge on the edge list of v then 

replace {v, w} as the current edge by the next edge on the list; 
(2b) else {{{v,w} is the last edge on the edge list of v )) begin 

make the first edge on the list the current one; 

perform cut{i) and change-value(i) for every child i of v; 

apply a price updating operation to v; 
end. 



Figure 2.10: The tree-push/update operation. 

implementation of the generic subroutine described in the previous section. How- 
ever, we replace the push/update operation with the tree-push/update operation 
described in Figure 2.10. 

A tree-push/update operation applies to a vertex v with a positive balance that 
is the root of a dynamic tree. There are two main cases. The first case occurs if the 
current edge {u, w} of v is eligible for a pushing operation. In this case we link these 
trees by making w the parent of v and do a send operation from v. The second 
case occurs if the edge {v, w} is not eligible for a pushing operation. In this case 
we update the current edge of v and update the price of v if necessary. If the price 
of v is increased, we cut all tree edges entering v, to maintain the invariant that all 
dynamic tree edges are admissible. 

It is important to realize that this algorithm stores values of the pseudoflow 
/ in two different ways. If {v,w} is an edge that is not a dynamic tree edge, 
f(v,w) is stored explicitly, with {v,w}. If {v,w} is a dynamic tree edge, with w 
the parent of v, then g(v) = u(v,w) — f(v,w) is stored implicitly in the dynamic 
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tree data structure. Whenever a tree edge (u, w) is cut, g(v) must be computed and 
/(u, w) updated to its current value. In addition, when the algorithm terminates, 
pseudoflow values must be computed for all edges remaining in dynamic trees. 

Two observations imply that the dynamic tree subroutine is correct. First, any 
edge {v,w} that is in a dynamic tree is admissible. By Lemma 2.6.5, in case (1) 
of tree-push/update, vertices v and w are in different trees, and the algorithm never 
attempts to link a dynamic tree to itself. Second, a vertex v that is not a tree 
root can have positive balance only in the middle of case (1) of a tree-push/update 
operation. To see this, note that only in this case does the algorithm add flow to a 
nonroot vertex, and this addition of flow is followed by a send operation that moves 
the nonroot excess of flow to one or more roots. 

Theorem 2.8.1 The dynamic tree implementation of Improve- Approximation sub- 
routine runs in 0(rcm log n) time. 

Proof: First we bound the number of link and cut operations. The number of 
link operations is at most Inm by a proof like that of Lemma 2.6.4. The number 
of cut operations is at most the number of link operations. The total number of 
link and cut operations is 0{nm). We charge 0(1) dynamic tree operations plus 
0(1) additional work to each cut and to each link operation. Since each dynamic 
tree operation requires 0(log n) work, the total amount of work charged to these 
operations is 0(nm\ogn). 

Consider a send operation. Work done at each iteration of the inner loop of send 
is charged to the cut operation performed at this iteration. If during an execution 
of the outer loop of send the inner loop is entered, the work done at this outer loop 
iteration (excluding the work done inside the inner loop which is already accounted 
for) is charged to the first cut operation performed by the inner loop at this iteration 
of the outer loop. If during an execution of the outer loop the inner loop is not 
entered, the send operation terminates after this iteration of the outer loop. We 
charge the work done by this iteration to the link operation performed in step (1) 
immediately before the send operation has been applied. This accounts for all the 
work done by send operations. 
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We account for the remaining work as follows. The work done at each execution 
of step (1) (except for the work done by send operations that is already accounted 
for) is charged to the link operation performed at this step. Each occurrence of 
step (2a) performs a constant amount of work, and the number of such occurrences 
is 0(nm) by a proof like that in Theorem 2.7.2. The work done by each pair of 
cut and change-value operations at step (2b) is charged to the cut operation. The 
work done by price updating operations adds up to 0(nm) by a proof like that 
in Theorem 2.7.2. The remaining work done at steps (2b) is 0(1) for each price 
updating operation, for a total of 0(n 2 ). 

The total amount of work charged to the cut and link operations is 0(nm log n). 
Other work amounts to 0(nm) for step (2a) and 0(nm) + 0(n 2 ) for step (2b). 
Running time of the subroutine is therefore 0(nm log n). | 

Corollary 2.8.2 An implementation of the minimum-cost flow algorithm that 
uses the dynamic tree version of the Improve- Approximation subroutine runs in 
(3(nmlog(n)log(nC)) time. 

Proof: Immediate from Theorems 2.4.1 and 2.8.1. I 



2.9 Blocking Improve- Approximation Subroutine 

In this section we present an alternative approach to the design of the Improve- 
Approximation subroutine. This approach is a generalization of the approach to 
the maximum flow problem due to Dinic and uses the concept of a blocking flow. 
The results presented in this section represent joint work with Robert Tarjan, who 
observed that the blocking flow approach can be extended to obtain an efficient 
implementations of the Improve- Approximation subroutine. 

A blocking flow is defined for flow networks (defined in Section 1.2), rather than 
for circulation networks. Recall that in flow networks all lower bounds on capacity 
are zero and there are two distinguished vertices, a source s and a sink t. A blocking 



2.9. BLOCKING IMPROVE-APPROXIMATION SUBROUTINE 71 

flow is a flow with no forward augmenting path, i.e. a flow for which for every path 
from s to t in the network contains a saturated edge. We use blocking flows only in 
the context of layered and acyclic networks. A network is acyclic if the underlying 
directed graph is acyclic. An equivalent definition is to say that a network is acyclic 
if its vertices can be labeled by integers in such a way that for every edge (u, w), we 
have label{v) > label(w). A network is layered [14] if its vertices can be labeled by 
integers in such a way that for every edge (v,w), we have label(v) — label(w) + 1. 

There are many algorithms for finding a blocking flow in a layered network. 
Although a layered network is a special case of an acyclic network, most of these 
algorithms work in the more general case as well, achieving the same complexity 
bounds. In particular, the algorithms described in [47,56,72] can be used to find a 
blocking flow in an acyclic network in 0(n 2 ) time, and the algorithm described in 
[67] can be used to find a blocking flow in 0{m log n) time. 

With a minor modification, the Shiloach-Vishkin algorithm [65] can be used 
to find blocking flows in acyclic networks in 0(n log n) PRAM time using 0(n) 
processors and 0{n 2 ) memory. A modification is necessary because, unlike a layered 
network, a general acyclic network does not have even/odd layers. One way to deal 
with this problem is to introduce a new vertex i( v<w ) for every edge (v, w), and replace 
the edge by two edges, (i>,i(„, w )) and (i( v<w ),w), both with the same capacity as the 
edge (v,w). The modified problem is clearly equivalent to the original problem, 
and in the modified network, the even/odd layers are well defined. The second way 
to deal with the problem is to modify the algorithm so that the flow pushed into 
a vertex during a parallel pulse is not processed until the end of the pulse. Both 
modifications are minor and do not affect the analysis. 

Galil's algorithm [27] finds blocking flows in layered networks in 0(n 2 / 3 m 2 / 3 ) 
time. This algorithm can also be modified to work on acyclic networks within the 
same time bound [28]. 

In this section, we define the notions of an admissible edge and an admissible 
graph in a slightly different way than in section 2.6. Given a pseudoflow / and 
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Procedure Improve- Approximation^ f , p, e); 

(( initialization)} 

V(u, w) G V X V do begin 

if c p (t;, w) > then f(v, w) *— l(v, w); 

if c p (v, w) < then /(u, u>) <— u(u, w); 
end; 

while / is not a circulation do block(f,p,€); 
return(/,p); 

end. 



Figure 2.11: The blocking flow Improve-Approximation subroutine. 

a price function p, we say that an edge (v,w) is admissible if (v,w) G Ef and 
c p (v,w) < 0. For a given price function p, the admissible graph Ga is the graph of 
all admissible edges: G A = (V, E A ), where E A = {(v,w) G Ef\c p (v,w) < 0}. These 
definitions are more natural in the context of blocking flows. As we have noticed at 
the end of Section 2.6, this definition can also be used in the previous sections as 
well. 

We show how to implement the Improve-Approximation subroutine using at 
most 3n blocking flow computations. The blocking flow variant of the subroutine, 
summarized in Figure 2.11, starts by bringing all edges in kilter. Then the subrou- 
tine repetively performs the block operation described in Figure 2.12 until / is a 
circulation. 

The block operation consists of two steps. In the first step, the vertices are 
partitioned into two sets, S and S. The set S contains the vertices reachable from 
some vertex with positive balance in the admissible graph G A , and the set S contains 
all remaining vertices. Then prices of vertices in S are increased by e/2. The second 
stage constructs an acyclic auxiliary network, finds a blocking flow /' in the auxiliary 
network, and augments the pseudoflow / by /' (i.e. f(v, w) <— f(v, w) + f'(v, w) for 
each (i>, w) G V x V). The auxiliary network G aux — (V \J{s, t}, E aux ) is constructed 
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Procedure block(f,p,e); 

« step 1)) 

S <— {wia-u such that b(v) > and w is reachable from v in Ga}] 

S^V-S; 

Vu e S do p(v) <- p(v) + e/2. 

(( step 2)) 

Construct auxiliary network G aux — (V\J{s,t},E aux ). 

Find a blocking flow /' in G aux . 

Augment /by /'. 

return(/, p); 

end. 



Figure 2.12: The block operation. The running time of the operation depends on the 
blocking flow subroutine used and on the details of implementation. 

from the admissible graph Ga by adding a source s, a sink t, and edges (s, v) of 
capacity bf(v) for every vertex v G V with bj{v) > and edges (w,t) of capacity 
—bf(w) for every vertex w G V with bj(w) < 0. The capacities of edges (v,w) G E A 
are defined to be rf(v,w). Note that augmenting / by a flow /' in the auxiliary 
network results in a valid pseudoflow. 

Correctness of the blocking flow subroutine follows from the following lemma. 

Lemma 2.9.1 At the beginning and at the end of each execution of block, pseud- 
oflow f is (e/ '2) -optimal with respect to p, S^ {v\bj(v) > 0}, and 5 2 {v\bf(v) < 
0}. At the beginning and at the end of steps (1) and (2), the admissible graph Ga 
is acyclic. 

Proof: We prove the lemma by induction on the number of block operations per- 
formed. The basis follows from the fact that after the initialization, we have Ea = 0. 
Now suppose that at the beginning of an execution of block the claim of the lemma 
holds. First we show that the (e/2)-optimality of / is preserved. Note that be- 
fore step (1) is executed, the reduced cost c p (v,w) > for every residual edge 
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(v,w) G S X S. Since the price change performed during step (1) decreases the 
reduced cost of such residual edges by e/2, the (e/2)-optimality constraints at these 
edges remain valid. Consider (v,w) E S x S such that (v,w) G G f . If (w,v) G G f , 
then the constraints for (w, v) remain valid by the previous case and the constraints 
for (v,w) remain valid by Lemma 2.3.3. If (w,v) <£ Gf, then f(v,w) = l(v,w) 
and the increase in the reduced cost c p (v,w) caused by the price change cannot 
invalidate the corresponding constraint. The reduced costs of the other residual 
edges do not change. The above observations imply that / is optimal after step (1). 
Step (2) changes the values of / on admissible edges only, and therefore / remains 
(e/2)-optimal after step (2). 

Suppose that at the beginning of step (1) the admissible graph is acyclic. The 
price increases performed at step (1) can introduce admissible edges going from S 
to S, but these price increases delete all admissible edges going from S to S. The 
admissible edges not going across the cut are not affected, and therefore Ga remains 
acyclic. Step (2) can delete admissible edges but cannot introduce new admissible 
edges, and therefore Ga remains acyclic after the step. 

Finally, suppose that after step (2) there is a path of admissible edges from a 
vertex with positive balance to a vertex with negative balance. Then there is a 
forward augmenting path with respect to /' in the auxiliary network, contradicting 
the fact that /' is a blocking flow. | 

The following lemma is a generalization of the lemma bounding the number of 
phases in a layered network algorithm for the maximum flow problem. This lemma 
is also similar to the lemma that bounds the number of maximum flow computations 
in a scaling step of a Bland- Jensen algorithm [8], and to Lemma 2.6.2. 

Lemma 2.9.2 The number of blocking flow computations in an execution of the 
blocking flow Improve- Approximation subroutine is at most Zn. 

Proof: Note that the balance of a vertex can become positive only during the 
initialization step of the subroutine and the subroutine terminates when there are 
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no vertices with a positive balance. Thus there must be a vertex v whose balance is 
positive during each execution of step (1) of phase operation. By Lemma 2.9.1, the 
price of v increases by e/2 during such a step, and does not change otherwise. By 
a similar argument, prices of vertices with a negative balance do not change. By 
an argument similar to the proof of Lemma 2.6.2, the price of v cannot increase by 
more then (3/2)e. Therefore the number of phases is at most 3rz. | 

Theorem 2.9.3 Let B(n,m) be the running time of an algorithm that finds a 
blocking flow in acyclic networks. Then the minimum-cost flow algorithm runs 
in 0(nB(n,m)\og(nC)) time. 

Proof: Immediate from Theorem 2.4.1 and Lemma 2.9.2. | 

Corollary 2.9.4 The minimum-cost flow problem can be solved in 
0(min(nm log(n) log(nC), n 5/3 m 2/3 log(nC), n 3 log(nC))) time. 

Proof: The lemma follows from Theorem 2.9.3 and on generalizations of blocking 
flow algorithms [47,27,67] to acyclic networks discussed above. | 



2.10 Parallel and Distributed Implementations 

In this section we discuss parallel and distributed implementations of the minimum- 
cost flow algorithms described earlier. Related work on parallel and distributed 
algorithms for network flow problems includes a paper of Shiloach and Vishkin [65], 
and a paper of Bertsekas [6]. Shiloach and Vishkin give an 0(n 2 log n) time parallel 
algorithm for the maximum flow problem, which can be combined with the cost 
scaling algorithms [62,8] to obtain an 0(n 3 log(n)(log C)) time parallel algorithm 
for the minimum-cost flow problem. Bertsekas [6] exhibits a "chaotic" algorithm 
for the minimum-cost flow problem that converges in a finite number of steps in a 
distributed model. 
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First we discuss two synchronous algorithms. One algorithm gives the better 
time bound, and the other gives the better memory bound. (In fact, if Conjecture 
1 is true, the time bounds for both algorithms are the same to within a constant 
factor, and the second algorithm is superior). The first algorithm, which we call the 
blocking algorithm, is obtained using the parallel version of the Shiloach-Vishkin 
blocking flow algorithm to implement the Improve- Approximation subroutine as 
described in Section 2.9. The second algorithm, which we call the FIFO algorithm, 
is obtained by using a parallel implementation of the first-in, first-out version of the 
Improve- Approximation subroutine. This parallel implementation can be obtained 
from the sequential implementation described in Section 2.7 in the same way as for 
the first-in, first-out maximum flow algorithm described in Section 1.4. 

First we state the performance of parallel versions of the algorithms. The bounds 
presented below hold in both PRAM [21] and DRAM [53] models of parallel com- 
putation. 

Theorem 2.10.1 The blocking algorithm runs in 0(ra 2 (logra)log(nC)) parallel time 
using 0{n) processors and 0{n 2 ) memory. 

Proof: Immediate from Theorem 2.4.1 and [65]. I 

Theorem 2.10.2 The FIFO algorithm runs in 0(n 3 (log n)(\ognC)) parallel time 
using 0(n) processors and 0{m) memory. 

Proof: Similar to the analysis of the parallel maximum flow algorithm (see Section 
1.6) using the 0(n 3 ) bound on the number of pulses of the Improve- Approximation 
subroutine instead if the 0(n 2 ) bound which holds in the maximum flow algorithm. 
I 

If Conjecture 1 holds, the running time bound for the FIFO algorithm improves 
by a factor of n, and the algorithm becomes superior to the blocking algorithm. 

The following two theorems give the performance of the two algorithms in the 
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synchronous distributed model of parallel computation [31,2]. In the statement of 
these theorems, A p denotes the degree of processor p in the network. 

Theorem 2.10.3 In the synchronous distributed model, the blocking algorithm runs 

in 0(n 2 (\ognC)) time using 0(nA p ) memory per processor p and 0(n 3 log(nC)) 
messages. 

Proof: Immediate from Theorem 2.4.1 and [65]. | 

Theorem 2.10.4 In the synchronous distributed model, the FIFO algorithm runs 

in 0(n 3 (\ognC)) time using 0(A P ) memory per processor p and 0(n 2 m\og(nC)) 
messages. 

Proof: Similar to the analysis of synchronous distributed version of the maximum 
flow algorithm. | 

If Conjecture 1 holds then the time complexity bound for the FIFO algorithm 
improves by a factor of n and the message complexity bound improves by a factor 
of m/n. 

Next we discuss implementations of the minimum-cost flow algorithm in the 
asynchronous distributed model of parallel computation [31,2]. One approach is to 
use the blocking or FIFO algorithm with the synchronization protocol of [2]. The 
resulting bounds on message and memory complexity of the algorithm are the same 
as in Theorems 2.10.3 and 2.10.4, and the resulting time bounds are greater then 
the bounds given by these theorems by a factor of log n. 

Another way to obtain an asynchronous implementation of the algorithm is to 
use an acknowledgment version of the first-in, first-out algorithm constructed as in 
the corresponding version of the maximum flow algorithm. The performance of the 
resulting algorithm is stated in the following theorem. If Conjecture 1 holds, the 
running time bound given by the theorem improves by a factor of n. 

Theorem 2.10.5 An asynchronous distributed version of the minimum-cost flow 
algorithm that uses the acknowledgment implementation of the Improve- Approxi- 
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motion subroutine runs in 0(n 3 log(nC)) time using 0(A P ) memory per processor 
p and 0(n 2 m\og{nC)) messages. 

Proof: Theorem 2.4.1, Lemma 2.7.3, and the proof of Theorem 1.6.3 are sufficient 
to prove the theorem if we can show how to detect termination of the Improve- 
Approximation subroutine. To detect termination, we use techniques from [31] and 
[2]; we assume familiarity with [31,2] in the following discussion. 

To detect the termination of the subroutine, we use the fact that the subroutine 
stops as soon as all vertices that had negative balance at the beginning of the exe- 
cution have nonnegative balance. At the beginning of the algorithm, we construct a 
spanning tree for the network. The root of the spanning tree plays a special role in 
the algorithm (without loss of generality, we assume that the network is connected). 
At the beginning of each execution of the Improve- Approximation subroutine, each 
processor records the sign of its initial balance, and sends the sign to the root of 
the tree. The root counts the number of vertices that had negative balance initially. 
When a balance of a vertex changes from negative to nonnegative, the vertex in- 
forms the root about this change, and the root updates its count of the number 
of vertices with negative balance. Consider the point when every processor has in- 
formed the root about the sign of its initial balance, and every vertex whose balance 
was negative initially has reported the change in sign. At this point, the root knows 
that the algorithm has terminated, so it broadcasts the termination massage. The 
number of messages related to termination detection is 0(n 2 ) for an execution of 
the subroutine, and the additional time is 0(n 2 ). The additional memory required 
is 0(1) for all vertices including the root. | 



2.11 Remarks 

The concluding remarks concern practicality and potential improvements and ex- 
tensions of the minimum- cost flow algorithm. 

Initial practical experience with the scaling algorithms for the minimum-cost 
flow problem turned out to be disappointing. As a result, the algorithms were con- 
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demned as impractical in the literature, and the simplex- based algorithms continued 
to be used to solve the minimum-cost flow problem in practice. This situation did 
not change even though the scaling algorithms were refined and new scaling al- 
gorithms were developed for the problem. Recently, however, Bland and Jensen 
conducted an extensive study [8] to compare an implementation of their scaling 
algorithm with simplex-based codes. They found their algorithm competitive, and 
concluded that the practicality of the scaling algorithms should be investigated fur- 
ther. These experimental results agree with the results of Bateson [3] who compared 
Gabow's scaling algorithm for general matching [26] with the Hungarian method 
[50]. 

We believe that our approach yields highly practical algorithms. Our belief 
is based on the work of Bland and Jensen mentioned above, as well as on the 
experience with implementations (both sequential and parallel) of the maximum 
flow algorithm described in Chapter 3 (also see [58]). The experience of Bertsekas 
and Tseng with an implementation of their algorithm [7] also supports our belief. 
Further experimental study is needed, however, before the practicality claim can be 
made with certainty. 

The approach presented in this chapter allows a great degree of flexibility. For 
example, the Improve- Approximation subroutine reduces the error parameter of 
approximation by a factor of 2. The subroutine can be modified to reduce the pa- 
rameter by a different factor. Although fine-tuning of this factor does not seem to 
improve the asymptotic running time of the algorithm, it does change the constant 
factors of the running time and therefore the practical performance of the algo- 
rithm. Also, as in the case of the maximum flow algorithm, a different ordering of 
the basic operations in the generic Improve- Approximation subroutine may result 
in a better performance. As we have mentioned in Section 2.7, we believe that 
the implementations of the Improve- Approximation subroutine that correspond to 
the first-in, first-out maximum flow algorithm and to the dynamic tree maximum 
flow algorithm achieve asymptotic running times equivalent to the corresponding 
maximum flow bounds. 

The algorithm can start with any initial price function such that the absolute 
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value of the difference between the prices of any pair of adjacent vertices is at 
most C (or even 0{n k C) for a constant k). In fact, in a practical implementation 
the initial price function should be obtained by using a shortest path subroutine. 
Improve- Approximation can also use the shortest path subroutine to update the 
price function at intermediate points of the execution. This use of the shortest 
path subroutine is similar to the use of breadth-first search in the maximum flow 
algorithm discussed in Section 1.7. Another heuristic improvement is to set e to the 
largest amount of violation of the complementary slackness conditions before each 
call to Improve- Approximation. 

Our algorithm works essentially by scaling costs. It would be interesting to see if 
a similar approach could be made to work by scaling capacities, or by scaling costs 
and capacities simultaneously. A possible approach is to order operations of the 
generic subroutine depending on the change in value (i.e. cost-capacity product) 
caused by these operations. 

Both the maximum flow and the shortest path problems are special cases of 
the minimum-cost flow problem. Our result implies that the complexity of the 
minimum-cost flow problem is not too far from the sum of the best upper bounds 
known for the two subproblems. An interesting extension would be to show how 
to solve the minimum-cost flow problem in 0(log(nC)) applications of any ("black- 
box") maximum flow and shortest path subroutines. Such a result would imply that 
the minimum-cost flow problem is not much harder then these two subproblems. 

Another interesting open question is the existence of a strongly polynomial 
0(n 3 log n) (or O (nm \og k n), if one is able to take advantage of sparsity) algo- 
rithm for the minimum-cost flow problem. A potential approach would combine 
the techniques of [70,59,23,30] with the techniques described in this chapter. 



Chapter 3 

Implementing Parallel Algorithms 



3.1 Introduction 

This chapter describes an implementation of the parallel maximum flow algorithm 
of section 1.6 on a Connection Machine 1 parallel supercomputer and presents ex- 
perimental results obtained using the implementation. The implementation makes 
an extensive use of the parallel prefix operations discussed in the next section. The 
use of parallel prefix operations was motivated by the work of Blelloch [9]. The 
material presented in this chapter represents joint research with Charles Leiserson. 
Most of the work, including all the experimental work, was done while the author 
was at Thinking Machines, Inc. 

In the previous two chapters we have discussed parallel and distributed ver- 
sions of the network flow algorithms. In this chapter we shall see how practical 
one of these algorithms is and how it compares with the underlying sequential algo- 
rithm. These experimental results are especially interesting because, until recently, 
large parallel machines were not available commercially, and practical experience 
with these machines is rare. Although important measures of parallel complexity 
(running time, processor requirements, communication, and memory) have been 
identified, the tradeoffs among these complexity measures are not yet understood 
completely, especially from the practical standpoint. As a result, analysis of a paral- 



1 Connection Machine is a trademark of Thinking Machines, Inc. 
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lei algorithm does not allow us to predict the practical performance of the algorithm 
as precisely as does the analysis of a sequential algorithm. 

Although implementation of parallel algorithms is an extensive area of research, 
the scarcity of practical experience with large scale parallel computers makes the 
choice among several possible implementations difficult. Another difficulty in pro- 
gramming parallel algorithms comes from incomplete understanding of the funda- 
mental parallel operations, i.e., operations used to express algorithms at a high 
level of abstraction that can be implemented efficiently on a wide class of parallel 
machines. 

We give an abstract model of the Connection Machine architecture that enables 
us to explain and justify the details of the implementation of the maximum flow 
algorithm. A detailed description of the design and the capabilities of the machine 
appears in the book of Hillis [41]. 

The Connection Machine can be viewed as a distributed memory parallel com- 
puter [53]. The machine consists of a large number of processors connected by a 
routing network. (The biggest machine available now has 64K processors, where 
K = 1024.) Each processor has a local memory. The local memory of a processor 
can be accessed by other processors via the routing network. The machine is a 
single-instruction, multiple-data machine [18]: the program is stored in a host com- 
puter which executes a sequential program containing parallel instructions. When a 
parallel instruction appears in the program, the host broadcasts it to all processors. 
Each processor, depending on its memory contents, either executes the instruction 
or remains idle. The operation of the machine is totally synchronous: the next 
instruction does not start until all processors have completed the execution of the 
current instruction. 

Due to its distributed nature, the memory organization of the machine is differ- 
ent from the shared memory organization used by some other parallel computers 
like the Ultracomputer [40]. In the Connection Machine, the interprocessor commu- 
nication is achieved by allowing a processor to access memory of any other processor 
using the network. Basic understanding of the memory system of the machine, in 
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particular the notions of a routing cycle and locality, is essential for the subsequent 
discussion. The locality of memory access is due to the fact that the time required 
by a processor to access its own memory is much less than the time required by a 
processor to access memory of another processor. More generally, the time required 
by a processor to access memory of a processor that is close (in the network) is 
greater than the time required to access memory of a processor that is far. This 
brings us to the notion of a routing cycle. Suppose that during the execution of an 
instruction several processors access memory of other processors. Since the opera- 
tion of the machine is synchronous, the longest memory access time determines the 
execution time of the instruction. A routing cycle is the amount of time it takes to 
execute such an instruction. It is important to realize that the routing cycle time 
varies depending on the interprocessor communication pattern and the size of data 
being accessed. 

When programming any computer at the user level, one wants to abstract the 
low-level details of the machine architecture. In particular, when programming a 
parallel computer, one wants to abstract the details of underlying routing network 
and routing algorithm. How is it possible to do so and still make use of locality? One 
answer is to use fundamental parallel operations like the parallel prefix operations 
discussed in the next section. 



3.2 Parallel Prefix Operations 

Parallel prefix operations have been recognized as fundamental parallel operations, 
and their implementation and use in describing parallel algorithms has been widely 
studied [9,11,21,49,51,53]. Given an associative binary operation "*" and a sequence 
of s numbers ai, a 2 , • • • « S5 the parallel prefix "*" operation maps this sequence into 
the sequence of s numbers a.\, a^ * a 2 , ■ ■ . , ct\ * a 2 * . . . * a s . 

We are especially interested in the following three operations. The first one 
is the prefix-add operation, obtained by using the plus function. This operation 
is illustrated in Figure 3.1a. The second operation is the prefix-copy operation, 
obtained by using the projection function Tt(x,y) = x. An application of this 
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(a) prefix-add[l 1111] = [12 3 4 5] 

(b) prefix-copy[l 2 3 4 5] = [11111] 

(c) prefix-min[4 5 2 13] = [44211] 

Figure 3.1: Parallel prefix operations. 

(a) seg-prefix-add([l 1 1][1 1][1][1 1 1]) = ([1 2 3][1 2][1][1 2 3]) 

(b) seg-prefix-copy([l 2 3] [4 5] [6] [7 8 9]) = ([1 1 1][4 4] [6] [7 7 7]) 

(c) 3eg-prefix-min([3 2 1}[4: 5}[Q][7 9 8]) = ([3 2 1][4 4] [6] [7 7 7]) 

Figure 3.2: Segmented parallel prefix operations. 

operation to a list of numbers results in copying the value of the first element of 
the list to all other elements of the list, as illustrated in Figure 3.1b. The third 
operation, prefix-min, is obtained using the minimum function. This operation is 
illustrated in Figure 3.1c. After this operation is applied, the last element of the 
list is equal to the smallest element of the input list. 

For our implementation we need an extension of parallel prefix operations. Given 
an associative binary operation "*", a segmented parallel prefix operation maps a 
list of sequences Si, 5 2 , . . . , Si into the list S[, S' 2 , . . . , 5/ where each sequence S- is 
obtained from the corresponding sequence 5, using the parallel prefix "*" opera- 
tion. The lengths of the input sequences 5 t - can be different. Segmented parallel 
prefix operations corresponding to addition, projection, and minimum functions are 
illustrated in Figure 3.2. 

Parallel suffix operations are defined similar to the parallel prefix operations, 
but enumerating sequences in reverse order. Segmented parallel suffix operations 
are illustrated in Figure 3.3. 

The parallel prefix operations, segmented or not, can be implemented efficiently 
on any distributed memory parallel machine such that a balanced binary tree can 
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(a) ^^^-^([l 1 1][1 1][1][1 1 1]) = ([3 2 1][2 1][1][3 2 1]) 

(b) seg-suffix-copy([l 2 3}[4 5][6][7 8 9]) = ([3 3 3] [5 5] [6] [9 9 9]) 

(c) seg-suffix-min([3 2 1][4 5] [6] [7 9 8]) = ([1 1 1][4 5] [6] [7 8 8]) 

Figure 3.3: Segmented parallel suffix operations. 

be embedded in the routing network of the machine. In particular, parallel prefix 
operations have been implemented on the Connection Machine [9]. To understand 
the use of these operations on the Connection Machine, one has to know that the 
processors on the machine are numbered. A sequence of numbers is represented 
in the machine by placing elements of the sequence in consequent processors. The 
processors containing the first and the last elements of the list are marked appro- 
priately. 

The efficiency of a parallel prefix computation depends of complexity of the 
binary operation "*" and on the size of data items. For simple operations like 
addition, projection, and minimum, the time required to perform parallel prefix 
operations is of the same order of magnitude as the time required for one routing 
cycle on data of the same size. Although the routing cycle time and time to perform 
a parallel prefix operation varies depending on the communication pattern and on 
the type of the parallel prefix operation, we shall refer to these times as to the 
routing cycle in our description of the implementation, since this gives a good 
enough intuitive measure of performance. The important point to remember is that 
the routing cycle is much greater then the local instruction execution time. 



3.3 Implementation Details 

We start the description of the implementation by describing the mapping of the 
input network into the machine. This mapping is due to Blelloch [9] who used it to 
implement several graph algorithms on the Connection Machine. 
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In our description of the mapping, it is important to realize that processors of 
the machine are linearly ordered in a way that allows efficient implementation of 
the segmented parallel prefix and suffix operations described in the previous section. 
Each vertex v is assigned a processor P v . Each undirected edge {v,w} is assigned 
two processors, P( V)U; ) and P( w>v y We call processors P( v , w ) and P( w , v ) pair processors. 
Vertex and edge processors are selected so that a vertex processor P v is immediately 
followed by the processors P( v<w ) that correspond to edges incident to v, and the 
last edge on the list is marked as such. This organization allows us to perform 
segmented parallel prefix and suffix operations on sequences consisting of vertices 
with their edges. It is important to realize that the information associated with 
an undirected edge {v,w} is stored in both P(„ )W ) and P( w<v ). In addition to the 
information associated with the edge {v , iu}, processor P(„ iU1 ) also stores the address 
of its pair processor P( WlV ), and the distance label values d(y) and d(w). 

The initialization step of the algorithm is easy to implement, so we shall concen- 
trate on implementation of the main loop (see section 1.6). The main loop consists 
of application of the pulse procedure until these are no active vertices. The pulse 
procedure is the only nontrivial part of the implementation. The implementation 
of this procedure is summarized on Figure 3.4. 

Steps (l)-(3) implement the first stage of the pulse procedure. Step 1 copies the 
value of the excess at each vertex processor P v , v E V — {.s,i}, to the edge processors 
P(v,w) using a seg-prefix-copy operation. The main purpose of step (2) is to distribute 
vertex excesses among outgoing edges. Our implementation favors edges at the end 
of edge lists. For each edge (v, w), we first compute the maximum amount that can 
be pushed through the edge. This amount is equal to Cf(v,w) if d(v) = d(w) + 1 
and to zero otherwise. Then, a seg-suffix-add operation is executed on these values. 
Now, each edge processor P( ViW ) has information about the excess to be pushed from 
v, the amount it can push, and the amount that can be pushed through the edges 
that follow (v ,w) on the edge list of v. This information is enough to compute the 
amount a(v,w) to be pushed through the edge (v,w). After the execution of the 
seg-suffix-add operation, each vertex processor P v contains the information about 



3.3. IMPLEMENTATION DETAILS 87 



Procedure Pulse. 

(1) For all v € V — {s,t} copy e(v) to all P( v<w ) using seg-prefix-copy operation 

(2) {( distribute excesses)) 

For all P( v>w ) compute, using seq-suffix-add, the amount that can be 

pushed to lower-labeled vertices through edges that follow (o, to) 
on the edge list of v. 

For all P( VtW ) compute the amount a(v, to) to be pushed from v to to. 

For all v € V compute the amount of excess that remains at v after the 
pushing. 

(3) (( push flow)) 

For all P(v,w) do if <r(v, to) > then begin 

/(», w) <- /(t>, to) + <r(», to); c/(», u>) ♦- c/(», to) - a{v, to); 

send a message containing o(v, to) to P( w ,v)i 
end. 
For all P( w , v ) that received message <r(t>, to) do begin 

/(to, v) *- /(to, v) - <r(o, to); c/(to, ©) ♦- c/(to, v) + o(v, to); 
end. 

(4) (( compute new distance labels)) 
For all P (ViW ) do 

if Cf(v, to) > then head-label^v, to) *- d(to) 
else head-label(v, to) <- 2n. 
For all v 6 V — {*, *} compute neto-d(») using seg-suffiz-min. 

(5) For all o € V - {«, *} copy new-a\v) to all P(o, to) using seg-prefix-copy. 

(6) ({ broadcast changed labels)) 

For all P( v , w ) such that v # {s, t} do 
if d(v) jL new - d(v) then 

set the value of d(v) at P( w , v ) to neto-d(»). 

(7) (( update excesses)) 

For all to e V do begin 

compute, using seg-suffix-add, the amount of flow neto-e(to) 

pushed into to; 

e(to) *- e(to) + neto — e(to); 
end. 
end. 



Figure 3.4: Implementation of the pulse procedure. 
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how much excess can be pushed from the vertex v at this pulse. The processor sets 
the value of its variable e(u) to the amount that will remain after the pushing. 

Finally, in step (3), all edge processors P( v , w ) for which the amount <r(v, w) com- 
puted in step (2) is positive, increase f(v,w) by a(v,w), decrease Cf(v,w) by the 
same amount, and send a message containing <r(v, w) to their dual processors P( w , v ) 
(i.e., writes to the agreed-upon location in the memory of P( w , v ))- Each proces- 
sor P( w , v ) that receives such a message decreases /(to, v) by a(v, w) and increases 
Cf(w, v) by the same amount. 

Steps (4)-(6) implement the second stage of the pulse procedure. First, each 
edge processor P( u , w ) sets its variable head-label to either the value of d(w) + 1 if the 
edge (v, w) is not saturated or to 2n if the edge is saturated. Recall that the value 
of d(w) is stored locally at P(u, to), so the above computation is performed without 
using the router of the machine. Also, note that the new label being computed 
is less than 2n, which explains the use of 2n as a value of head-label for saturated 
edges. Next, the seg-suffix-min operation is performed on the head-label variable, 
and as a result each vertex processor P„ contains the new value of d(v). In step (5), 
all vertex processors except P, and P t copy this new value to their edge processors 
using a seg-prefix-eopy operation. In step (6), each edge processor P( v , w ) checks if the 
new value of d(v) is different from the old value and if it is, the processor updates 
its value of d(v) and the value of d(v) of its pair processor P( w , v ) (via the routing 
network). 

Step (7) implements the third stage of the pulse subroutine. The step performs 
the seg-suffix-add operation to compute, for each vertex v, the amount of flow pushed 
into v during step (3). This amount is added to the excess e(v) updated at step (2) 
to account for the flow pushed out of v. 

The running time of each step is slightly greater than a routing cycle of the 
machine, because the running time of each step is dominated by a segmented prefix 
or suffix operation (steps (1), (2), (4), (5), (7)) or by a message routing operation 
(steps (3), (6)). The overall running time of the pulse procedure is roughly seven 
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routing cycles of the machine, which results in a tight inner loop: the key factor in 
the performance of the implementation. 

We conclude this section with a discussion of the implementation. This imple- 
mentation is done in the spirit of data parallelism [42]: each item of data is assigned 
a processor which operates on this item. This approach is natural for the Connec- 
tion Machine with its fine-grained parallelism, characterized by a low memory to 
processor ratio. 

What are the tradeoffs of the implementation? The advantages of the implemen- 
tation are a tight inner loop and a simple but effective use of locality through parallel 
prefix operations. The disadvantage is a relatively low processor utilization. As we 
have mentioned in Section 1.6, the algorithm can be implemented with 0{n) pro- 
cessors instead of 0(m) processors used by the implementation we have described. 
The alternative implementation achieves better processor utilization. This imple- 
mentation, however, requires processor scheduling, which introduces extra overhead 
and makes the use of locality difficult. 

The above implementation is designed for the Connection Machine architecture, 
but it is also good for other kinds of architectures. For example, fetch-and-add 
operations can be used to implement parallel prefix operations on the Ultracomputer 

[40]. 



3.4 Experimental Results 

We start this section by describing the actual implementation, including the heuris- 
tics used to speed up the algorithm. Then we describe the generator of examples 
used to obtain the input data. Finally, we present and discuss the experimental 
results. 

The actual code implemented the minimum-cut variation of the maximum flow 
algorithm, as described in Section 1.7 and in [37]. The main difference of the mincut 
algorithm is that a vertex becomes "dead" as soon as its label reaches or exceeds n 
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(in this case, the sink is not reachable from the vertex). The implementation uses 
two heuristic improvements. The first improvement is to declare a vertex "dead" 
as soon as the value of its label exceeds n minus the number of "dead" vertices. 
The second improvement is to do breadth-first search backwards from the sink in 
the residual graph to update the distance labels to be the exact distances to the 
sink. The frequency of the breadth-first search updates is determined by the size of 
the graph and the number of relabeling operations since the previous breadth-first 
search update. 

The experiments have two goals: determine the practicality of the parallel im- 
plementation and compare the parallel implementation with the sequential imple- 
mentation of the same algorithm. In addition to the parallel implementation on the 
32K Connection Machine, the corresponding first-in, first-out algorithm was imple- 
mented on the Symbolics 3600 series Lisp Machine. The parallel implementation 
is written in *LISP and the sequential implementation in LISP. Both implemen- 
tations use the heuristics described above. The only difference is in the tuning 
of the breadth-first search update frequency. Different tuning is needed because 
of the different relative costs of breadth-first search in the parallel and sequential 
implementations . 

Our experiments are conducted on large and difficult examples of the minimum- 
cut problem. Generating such examples is a non-trivial task. Several methods we 
tried either generated trivial examples (where either the source or the sink is on one 
side of the cut and all other vertices are on the other side) or simple examples. On 
these examples, the programs run extremely fast, and the distribution of operation 
performed during an execution is different from what the analysis predicts: relabel- 
ing operations dominate the sequential running time. The method described below 
produces examples that take longer time and produce the distribution of operations 
that agrees with the worst-case analysis, with nonsaturating pushes dominating the 
sequential running time. 

To see how examples are generated, imagine an infinite pipe with a mesh drawn 
on it. Suppose the pipe goes from west to east. The distance from a vertex of the 
mesh to the nearest neighbor in both the horizontal and vertical directions is one, 
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and circumference of the pipe is D. First we construct a graph on the vertices of 
the mesh. In this graph, every vertex will have degree 4A, were A is assumed to be 
less then 25/2. To construct the graph, we connect each vertex v by a directed edge 
to all vertices u> within distance A due east, west, south, and north from v. The 
capacity of an edge (t>, w) is denned depending on the distance x between v and 
w. The capacity is selected at random from a uniform distribution on the interval 
[0,2^). 

To complete the construction, we introduce a source s and a sink t. Then we 
"cut" the pipe by two planes perpendicular to the axis of the pipe. The cutting 
planes are distance L apart. We disregard the portion of the pipe that is not between 
the planes. The edges cut by one plane are connected to the source preserving their 
directions and capacities; the edges cut by the other plane are connected to the 
sink in the similar manner. Note that since a single vertex may have several of its 
edges cut by a plane, a single vertex may have several edges connecting it with the 
source or the sink. These multiple edges are merged together, and the capacity of 
the resulting edge is set to the sum of capacities of the merged edges. 

Why is an average problem constructed by such a generator hard? In the net- 
work, short source-to-sink paths have small residual capacity, and long source-to- 
sink paths have large residual capacity. The algorithm tends to processes short 
paths first. Saturating short paths does not, however, saturate the cut, so the 
algorithm must continue. 

In our experiments, the values of parameters are selected using the following 
formulas: L = D = Wn\ , A = [(y/n- 1)/2J. The value of A determines the graph 
density and is selected to produce networks of moderate density. The values of n 
are selected so that y/n is integer and the size of the problem, i.e., the number of 
processors required by our implementation, is close to a multiple of the machine 
size. Two sequences of examples were evaluated. Tables 3.1a and 3.1b summarize 
the experimental data. 

As one can see from the table, the size of some examples is greater than 32K, the 
size of the Connection Machine used in our experiments. In these cases, the built-in 
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(a) 







Running Time 




n 


Sise 


LUp Machine 


Connection Machine 


Speedup 


531 


15,435 


56.8 


1.98 


29 


843 


31,699 


195.3 


3.24 


60 


1371 


65,307 


389.8 


4.42 


88 


2211 


130,803 


1111.3 


9.98 


111 



(b) 







Running Time 




n 


Size 


Liap Machine 


Connection Machine 


Speedup 


531 


15,435 


81.2 


2.33 


35 


843 


31,699 


142.7 


2.93 


49 


1371 


65,307 


303.1 


4.47 


68 


2211 


130,803 


1734.5 


9.58 


181 



Table 3.1: Experimental results for the first (a) and the second (b) sequence of exam- 
ples. All running times are in seconds and exclude paging time. 



virtual processor mechanism of the Connection Machine is used. This mechanism 
is simple but effective, and dpes not introduce any overhead. Our biggest examples 
require each real processor to serve four virtual processors. The first example in each 
sequence is much smaller than the machine size of 32K, so the processor utilization 
is low and the speedup compared to the sequential implementation is low as well. 

The running times given in the table are in seconds and do not include paging 
time, which is a problem for the Lisp Machine on larger examples. The running 
times for the parallel implementation are quite good: under ten seconds even for 
the largest examples. 

How does the parallel implementation compare with the sequential implementa- 
tion? The speedup varies from under 30 on small examples (when not all processors 
are utilized) to well over a hundred on the largest examples. In general, the speedup 
increases with the problem size. Two factors are responsible for this increase. First, 
bigger problems exhibit larger amounts of parallelism. Second, a larger number 
of virtual processors per real processor increases the locality of data and therefore 
decreases the relative cost of parallel prefix operations. 
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Several simple enhancements in Connection Machine systems software will sig- 
nificantly improve the performance of the parallel implementation. With these 
enhancements, we expect a 600-700 times speedup for the larger 64K Connection 
Machine on networks of size around 512K. 

In conclusion, we would like to emphasize that our experimental results are 
highly sensitive to the conditions of the experiment. For example, using higher- 
degree graphs would probably increase the speedup achieved by the parallel im- 
plementation, and using lower-degree graphs would decrease the speedup. In the 
low-degree case, the Ahuja-Orlin sequential algorithm [1] mentioned in Section 1.1 
should produce better results then the sequential algorithm that we have imple- 
mented. Since many maximum flow problems that appear in practice are much 
smaller and much simpler than the examples we have generated, most network flow 
algorithms produce satisfactory results. In Chapter 2 we have seen, however, that 
minimum-cost flow problems can be solved by iterating a generalized version of the 
maximum flow algorithm. Since large and hard instances of the minimum-cost flow 
problem do appear in practice, fast maximum flow algorithms are very important 
in this context. 
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Chapter 4 

Parallel Symmetry-Breaking 



4.1 Introduction 



In the previous chapters we have seen the need for new methods in design and 
implementation of efficient parallel algorithms. One way to identify these methods 
is to study problems for which simple but efficient sequential algorithms exist, but 
which appear to be much harder to solve efficiently in a parallel framework. A 
known example of a problem with a trivial sequential algorithm which is hard to 
solve in parallel is the problem of finding a maximal independent set in a graph 
[73]. This problem was shown to be in the class NC of problems which can be 
solved in polylogarithmic time using a polynomial number of processors by Karp 
and Wigderson [46]. A simple randomized algorithm for the problem is due to Luby 
[55]. Recently M. Goldberg and Spencer [39] gave a deterministic algorithm for the 
problem that runs in polylogarithmic time using a linear number of processors. 

The study of the maximal independent set problem shows the importance of 
techniques for breaking symmetry in parallel. The symmetry-breaking comes up 
in many other parallel algorithms as well. In many cases, however, it is enough 
to be able to break symmetry in special kinds of graphs. The performance of the 
resulting algorithm improves if we can solve the special case of symmetry-breaking 
more efficiently. 

95 
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In this chapter we present a technique for breaking symmetry. In particular, 
we give an 0(lg*n) time algorithm to 3-color a rooted tree. Our technique can be 
viewed as a generalization of the deterministic coin-flipping technique of Cole and 
Vishkin [13]. To show the usefulness of our technique, we present the following 
algorithms. All of the algorithms presented use a linear number of processors. 

• For graphs whose maximum degree is a constant A, we give an 0( A log A lg*n) 
algorithm for (A+l)-coloring and for finding a maximal independent set on 
an EREW PRAM. 

• We give an algorithm to 7-color a planar graph. Both this algorithm and the 
maximal independent set (for planar graphs) algorithm based on it run in 
0(lognlg*n) time on a CRCW PRAM and in O(log 2 n) time on an EREW 
PRAM. We also give an 0(log 3 n lg*n) CRCW algorithm to 5-color a planar 
graph. 

• We give an 0(log n lg*n) time algorithm for finding a maximal matching in a 
planar graph on a CRCW PRAM. 

The results stated above improve the running time and processor bounds for the 
respective problems. The fastest previously known algorithm for (A + l)-coloring 
[55], in the case of constant-degree graphs, runs in O(logn) time, and the deter- 
ministic version of this algorithm requires n 3 processors. The 5-coloring algorithm 
for planar graphs, due to Boyar and Karloff, [10] runs in 0(log 3 n) time, but the 
deterministic version of this algorithm requires n 3 processors. The 0(log 3 n) run- 
ning time of the maximal matching algorithm due to Israeli and Shiloach [44] can 
be reduced to 0(log 2 n) in the restricted case of planar graphs, but our algorithm 
is faster. 

Although in this chapter we have limited ourselves to the application of our 
techniques to the PRAM model of computation, the same techniques apply to a 
distributed model of computation [31,2]. Moreover, the 0(lg*n) lower bound given 
by Linial [54] for the maximal independent set problem on a chain in the distributed 
model shows that our symmetry-breaking technique is optimal in this model. 
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The fact that a rooted tree can be 3-colored in 0(lg*n) time raises the question 
whether a rooted tree can be 2-colored within the same time complexity. We answer 
this question by giving an ft(log n/ loglog n) lower bound for 2-coloring a rooted tree. 
We also present an ft(log n/ loglog n) lower bound for finding a maximal independent 
set in a general graph, thus answering the question posed by Luby [55]. 

This chapter represents joint research with Plotkin [34,35]. Results similar to 
many of the results presented in this chapter were obtained independently by Shan- 
non [63]. A joint paper will appear in [36]. 

4.2 Definitions and Notation 

This section describes our assumptions about the computational model and intro- 
duces the notation used throughout the chapter. As in the earlier chapters, we use 
n to denote the number of vertices and m to denote the number of edges in a graph. 
We use A to denote the maximum degree of the graph. 

Given a graph G = ( V, E), we say that a subset of vertices I € V is independent 
if no two vertices in I are adjacent. A coloring of a graph G is an assignment 
C : V — ► N of positive integers (colors) to the vertices of the graph. A coloring is 
valid if no two adjacent vertices have the same color. The I th bit in the color of a 
vertex v is denoted by C v (i). A subset of edges M € E is a matching if any two 
distinct edges in M have no vertices in common. 

The following problems are discussed in this chapter: 

• The vertex-coloring (VC) problem: find a valid coloring of a given graph that 
uses at most A+l colors. 

• The maximal independent set (MIS) problem: find a maximal independent 
set of vertices in a given graph. 

• The maximal matching (MM) problem: find a maximal matching in a given 
graph. 
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We make a distinction between unrooted and rooted trees. In a rooted tree, each 
nonroot vertex knows which of its neighbors is its parent. 

The following notation is used: 

lga: = log 2 x 

lgW x = lg x 

l g (0:r = lglg^-^x 

lg*x = min { , w . )a! < 2} (0 

We assume a PRAM model of computation where each processor is capable of ex- 
ecuting simple word and bit operations. The word width is assumed to be O(log n). 
The word operations we use include bit-wise boolean operations, integer compar- 
isons, and unary- to-binary conversion. In addition, we assume that each processor 
has a unique identification number O(logn) bits wide, which we denote by PE-ID. 
We use the exclusive-read, exclusive-write (EREW) PRAM, the concurrent-read, 
exclusive-write (CREW) PRAM, or the concurrent-read, concurrent-write (CRCW) 
PRAM as appropriate. For a discussion of the PRAM model of computation, see 
[21]. All lower bounds are proved for a CRCW PRAM with a polynomial number 
of processors. 



4.3 Coloring Rooted Trees 

This section describes an 0(lg*n) time algorithm for 3-coloring rooted trees. First 
we describe an 0(lg*n) time algorithm for 6-coloring rooted trees. Then we show 
how to transform a 6-coloring of a rooted tree into a 3-coloring in constant time. 

The procedure 6- Color- Rooted- Tree is shown in Figure 4.1. This procedure ac- 
cepts a rooted tree T = (V, E) and 6-colors it in time 0(lg*n). Starting from the 
valid coloring given by the processor ID's, the procedure iteratively reduces the 
number of bits in the color descriptions by recoloring each nonroot vertex v with 
the color obtained by concatenating the index of a bit in which C„ differs from 
Cjather(v) and the value of this bit. The root r appends C,[0] to 0. 
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Procedure 6-Color-Rooted-Tree. 






£<-flgn|; 






for all v € V in parallel do C v <- 


PE-ID(t>); 


(( initial coloring)) 


while L > flgX + 1"| for all » G V 


in parallel do begin 


if v is the root then begin 






**<-<); 






K <- C„(0); 






end; 






else begin 






*„ <- min^ | c v (i)^c fath9Tiv) ( 


••)}(•); 




K *- C v (i v ); 






end; 






C v *— b v t v ', 






end; 






end. 







Figure 4.1: The coloring algorithm for rooted trees 



Theorem 4.3.1 The algorithm 6- Color- Rooted- Tree produces a valid 6- coloring of 
a tree in 0(lg*n) time on a CREW PRAM using 0(n) processors. 

Proof: First we prove by induction that the coloring computed by the algorithm is 
valid, and then we prove the upper bound on the execution time. 

Assume that the coloring C is valid at the beginning of an iteration. We shall 
show that the coloring at the end of the iteration is also valid. Let t; and w be two 
adjacent vertices; without loss of generality assume that t; is the father of w. By 
the algorithm, w chooses some index t such that C v (i) ^ C w (i) and t; chooses some 
index j such that C v (j) £ C father ( v )(j). The new color of w is (i,C w (i)) and the 
new color of v is (J, C v (j)). If * ^ j, the new colors are different and we are done. 
On the other hand, if t = j, then C v (i) ^ C„,(t) and again the colors are different. 
Hence, the validity of the coloring is preserved. 

Now we show that the algorithm terminates after 0(lg*n) iterations. Let Lk 
denote the number of bits in the representation of colors after Jb iterations. For 
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k = 1 we have 

L x = flgXl+1 
< 2[lgLl 

if rig l] >i. 

Assume for some k we have .£*_! < 2flg (fc-1) L] and (lg (fc) £] > 2. Then 

L* = Pg2W|+l 

< [lg(21g( fe - 1 )L)l+l 

< 2 fig W I] 

Therefore, as long as flg^ L] > 2, 

£*<2flgWLl. 

Hence, the number of bits in the representation of colors Lk decreases until, after 
0(lg*n) iterations, flg (At) L] becomes 1 and L k reaches the value of 3 (the solution 
of L = flglr] + 1). Another iteration of the algorithm produces a 6-coloring: 3 
possible values of the index t'„ and 2 possible values of the bit b v . The algorithm 
terminates at this point. 

We use the concurrent-read capability to enable all the sons of a vertex v to 
access the color of t; simultaneously; no concurrent-write capabilities are required. 
For constant-degree trees the concurrent-read capability is not needed. | 

As we have shown, a rooted tree can be 6-colored quickly. A natural question to 
ask at this point is whether one can use fewer colors and still stay within the same 
complexity bounds. The following theorem answers this question. 

Theorem 4.3.2 A rooted tree can be S- colored in 0(lg*n) CREW PRAM time 
using 0(n) processors. 

Proof: The algorithm S- Color-Rooted- Tree presented in Figure 4.2 starts by using 
the previously described algorithm to 6-color the tree and then recolors it in 3 colors 
in constant time. 
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The algorithm recolors the vertices colored with bad colors 3, 4, and 5, into good 
colors 0, 1, 2 as follows. First, each vertex is recolored in the color of its father, so 
that any two vertices with the same father have the same color. The root, which 
has no father, recolors itself with a color different from its current color. Next, 
the algorithm removes the color from every vertex that has a bad color and has a 
neighbor with a good color. These vertices become uncolored. Every vertex v that 
still has a color C v is recolored in the color C v mod 3; this gets rid of the remaining 
bad colors. Note that this coloring has the property that for any vertex u, all of the 
sons of v that are colored must have the same color. 

The resulting coloring is valid, but not all vertices are colored. By the construc- 
tion, every uncolored vertex has at least one colored neighbor. Therefore, if there 
are two vertices v and w, such that v = fathtriw) and both vertices are uncolored, 
then father(v) is colored and aons(w) are colored too. The algorithm colors v with 
a color different from C ton ^ v ) and from Cf at her(v)- Such a color always exists because 
there are 3 different colors to choose from and all the colored sons of v have the 
same color. Finally, the algorithm colors to with a color different from both C„ 
and Cgoniiw). Every step of the 8-Color-Rooted-Tree algorithm can be executed in 
constant time except for the first step, in which we color the tree with 6 colors. 
Hence, the total running time of the algorithm is 0(lg*n). | 

Remark: Theorem 4.3.2 holds for a more general class of graphs containing all 
directed graphs such that out-degree of every vertex is at most one. The same 
proof applies with one small modification. Since the concept of a root makes no 
sense for graphs which are not trees, the proof should use a more general concept 
of a vertex with zero out-degree instead. 

Any tree can be 2-colored. In fact, it is easy to 2-color a tree in polylogarithmic 
time. For example, one can use treefix operations [53] to compute the distance from 
each vertex to the root, and color even level vertices with one color and odd level 
vertices with the other color. It is harder to find a 2-coloring of a rooted tree in 
parallel, however, than it is to find a 3-coloring of a rooted tree. In section 4.6 
we show a lower bound of fl(log n/ loglog n) on 2-coloring of a directed list by a 
CRCW PRAM with a polynomial number of processors, which implies the same 
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Procedure 8-Color-Rooted-Tree. 
C <- 6-Color-Rooted-Tree (V,E); 

for all t; e V, v / root in parallel do C v «- C/ 0<fcer ( v ); 
Croot «- min.-6{o,i,2}-{c. onj(POOt) } (*); 
*i «- {t> | C v < 2}; 
F 2 <-F-Vi; 

V <- {v | w e V 2 and 3(v,w) G £,«; € Vi}; 
for all v g V - V in parallel do C v <- C v mod 3; 
for all v € V in parallel do C v «- uncolored; 
for all t; 6 V in parallel do 
if father(v) £ V then begin 

C„ «- mii^ci.aMc^.jMCy.rt^.,} Oi 
V - F' - {»}; 
end; 
for all v € V in parallel do 

C v *- min,€ { o, 1 ,2 } _{ C . OBK . ) }-{c /a , her(B) }(t); 
end. 



Figure 4.2: The 3-coloring algorithm for rooted trees 
lower bound for rooted trees. 

4.4 Coloring Constant-Degree Graphs 

The method for coloring rooted trees described in the previous section is a gener- 
alization of the deterministic coin-flipping technique described in [13]. The method 
can be generalized even further [35] to color constant-degree graphs in a constant 
number of colors. In the generalized algorithm, the current color of a vertex is 
replaced by a new color obtained by looking at each neighbor, appending the index 
of a bit in which the current color of the vertex is different from the neighbor's 
color to the value of the bit in the vertex color, and concatenating the resulting 
strings. This algorithm runs in O(lg'n) time, but the number of colors, although 
constant, is exponential in the degree of the graph. In this section we show how to 
use the procedure S- Color- Rooted- Tree described in the previous section to color a 
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Procedure Find-Forest(V,E). 

iZ<-0; 

for all v £ V in parallel do {{ construct the forest - the first step}) 
if PE-ID(v) is not a local maximum then begin 
e v *- (v,w) s.t. (v,w) € E and PE-ID(w) = max{PE-ID(«)|(», «) € E}; 
E' <-E' U e„; 
end; 

else R <- R U v 
for all v € R in parallel do (( get rid of zero-depth trees - the second step)) 
if j!(w, w) € E' and 3(v, «/) G E then 
£'<-£'U (»,«/); 
return (E'); ({ the edges of the forest)) 
end. 



Figure 4.3: The spanning forest algorithm 

constant-degree graph with (A+l) colors, where A is the maximum degree of the 
graph, in 0( A 2 log( A) lg*n) time. 

First, we describe how to find in constant time a forest in a given graph such that 
each vertex with nonzero degree in the graph has nonzero degree in the forest. The 
removal of the edges of the forest decreases the maximum degree of the remaining 
graph (unless the maximum degree of the graph is zero). We shall use this property 
(later) to decompose the edges into A sets, each set inducing a forest on the vertices 
of the graph. The procedure Find-Forest (see Figure 4.3) constructs such a forest. 

The procedure has two steps. In the first step each vertex compares the ID's of 
its neighbors with its own ID. A vertex that does not have the maximum processor 
ID among its neighbors chooses an edge that connects it to the neighbor with the 
largest processor ID. The graph induced by the chosen edges is a forest (this graph 
has no cycles) and the vertices with the highest processor ID's among their neighbors 
- local maxima - are roots of the forest. In the second step each root with no sons 
chooses an edge that connects it to one of its neighbors. The roots are local maxima 
and are therefore independent. Hence, no new cycles are introduced into the graph 
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Procedure Color-Constant-Degree-Grttph. 
E' *-E; 
t<-0; 

while £' ^ do begin (( the first phase)) 
Ei <- Find-Fore$Hy,E r ); 
E' <-E'- Ec, 
i «- i + 1; 
end; 
for all v EV in parallel do (( initial coloring)) 

C(v) +- 1; 
for t <- i - 1 to do begin ({ the second phase)) 
C *- 3-Color-Rooted Tree (V, E t ); 
E'*-E' + E { ; 
for k *- 1 to 3 do 

for j *- 1 to A + l do begin 
V *- V; 
for all v € V in parallel do 

if C(v) m j and C'(v) = k then begin 

C(v) <- maxizu lt 2,..AH}-{c(u,) \ («,«)€£'}}(*); 
V - V - {«}; 
end; 
end; 
end; 
end. 



Figure 4.4: The algorithm for coloring constant-degree graphs 

induced by the chosen edges. 

The algorithm Color- Constant- Degree- Graph that colors constant-degree graph 
with (A+l) colors is presented in Figure 4.4. The algorithm consists of two phases. 
In the first phase we iteratively call the Find-Forest procedure, each time removing 
the edges of the constructed forest. This phase continues until no edges remain, at 
which point we color all the vertices with one color. 

In the second phase we iteratively add the edges of the forests back into the 
graph, each time recoloring the vertices to maintain a consistent coloring. At the 
beginning of each iteration of this phase, the edges of the current forest (E') are 
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added, making the existing (A + l)-coloring inconsistent. This forest is colored with 
3 colors using the S- Color- Rooted-Tree procedure. Now, each vertex has two colors 
- one from the coloring at the previous iteration and one from the coloring of the 
forest. The pairs of colors form a valid 3(A+l)-coloring of the graph. The iteration 
finishes by enumerating the color classes, recoloring each vertex of the current color 
with a color from {0, . . . , A} that is different from the colors of its neighbors (note 
that we can recolor all the vertices of the same color in parallel because they are 
independent). 

In the proof of the following theorem, we assume that A = O(log n). Under this 
assumption, we can use word operations to implement, in constant time, union and 
membership testing operations on sets of size A + 1. 

Theorem 4.4.1 The algorithm Color- Constant-Degree- Graph runs in 
0( A log A( A + lg*n)) time and colors the graph with (A+l) colors. 

Proof: At each iteration all edges of a spanning forest are removed. From the above 
discussion it follows that each vertex that still has neighbors in the beginning of 
an iteration has at least one edge removed during that iteration and therefore its 
degree decreases. Hence, the first phase of the algorithm terminates in at most A 
iterations. 

The second phase terminates in at most A iterations as well. Each iteration 
consists of two stages. First, the current forest is colored using procedure S- Color- 
Rooted- Tree, which takes, by theorem 4.3.2, O(log A lg*n) time on an EREW PRAM 
(the log A factor appears because we do not use the concurrent-read capability). 
Now we iterate over all the colors. Each iteration can be done in O(logA) time 
using word operations. Hence, one iteration of the second phase takes 0(log A lg*n+ 
A log A) time, leading to an overall 0( A log A( A+lg*n)) running time on an EREW 
PRAM. | 

Having a (A+l)-coloring of a graph enables us to find an MIS in this graph. The 
following theorem states this fact formally. (We refer to the algorithm described in 
the proof as Constant- Degree- MIS in the subsequent sections.) 
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Theorem 4.4.2 An MIS in constant- degree graphs can be found in 0(lg*n) time 
on an EREW PRAM using 0(n) processors. 

Proof: After coloring the graph in a constant number of colors using the procedure 
Color- Constant-Degree- Graph, one can find an MIS by iterating over the colors, tak- 
ing all the remaining vertices of the current color, adding them to the independent 
set, and removing them and all their neighbors from the graph. By theorem 4.4.1, 
the coloring of a constant-degree graph takes 0(lg*n) time on an EREW PRAM. 
The selection of all vertices with a specific color and the removal of all neighbors of 
the selected vertices takes constant time. | 

The proofs of theorems 4.4.1 and 4.4.2 also imply that the algorithms Color- 
Constant- Degree- Graph and Constant- Degree- MIS have a polylogarithmic running 
times for graphs with polylogarithmic maximum degrees. In this case, however, the 
assumption that the word size is greater than A is unreasonable, so the running 
time of the algorithms becomes 0(A(A 2 + log Alg*n)). 

The above algorithms can be implemented in the distributed model of com- 
putation [31,2], where processors have fixed connections determined by the input 
graph. The algorithms in the distributed model achieve the same 0(lg*n) bound 
as in the EREW PRAM model. Linial has recently shown [54] that fi(lg*n) time 
is required in the distributed model to find a maximal independent set on a chain. 
Our algorithms are therefore optimal (to within a constant factor) in the distributed 
model. 

4.5 Algorithms for Planar Graphs 

In this section we study the problem of coloring a planar graph and the problem of 
finding a maximal matching in a planar graph. 

Any planar graph can be 4-colored. Linear-time sequential algorithms, however, 
are known only for 5-coloring planar graphs. In this section we describe a simple and 
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Procedure 7-Color-Planar-Graph. 

i *- 0; 

(( first stage)) 

while V / for all t? € V in parallel do begin 
if Degree(t;) < 6 then begin 
Vi<-Vi + {v); 
V't-V- {»}; 
end; 

i *~ i + 1; 
end; 

(( second stage )) 
for t <— t - 1 to do 

while V; ?£ do begin 

Ei «- {(«,«>) | v, w eVii (v,w)€E); 
I*- Constant-Degme-MIS(yi,Ei); 
for all v € J in parallel do 

C v <- max,- 6 {{ 1 .„7}_.{ C . |l , | u, € v';(«, w )€K} }(*); 
V'*-V' + /; 
K^-Vi-/; 
end; 
end. 



Figure 4.5: The algorithm for 7-coloring of planar graphs 

efficient parallel algorithm that 7-colors a planar graph, and show how to construct 
a more complicated parallel algorithm to 5-color a planar graph. 

First we describe an algorithm for 7-coloring of planar graphs. The algorithm, 
called 7- Color- Planar- Graph, is shown in Figure 4.5. The algorithm consists of two 
stages. In the first stage, we iteratively partition the vertices of the graph into 
layers. At each iteration we create a new layer consisting of all vertices of the graph 
with degree 6 or less and delete these vertices from the graph. 

The second stage returns the layers to the graph in the order opposite to the 
order in which the layers are removed. After a layer is returned, it is 7-colored in a 
way consistent with the coloring of the layers which have been returned and colored 
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in the previous iterations. Note that all the vertices of the returned layer have a 
degree of at most 6 in the current graph. 

The layer is colored by iteratively applying the Constant- Degree- MIS procedure 
to find an MIS in the subgraph induced by the uncolored vertices of the layer, and 
coloring each of the selected vertices in a color different from its colored neighbors. 
Since the uncolored vertices have a degree of at most 6 in the current graph, we 
never need more than 7 colors. 

Theorem 4.5.1 The algorithm 7- Color- Planar- Graph runs in O(lognlg*n) time 
on a CRCW PRAM and in O(log 2 n) time on an EREW PRAM. 

Proof: In a planar graph, at least a constant fraction (1/7) of the vertices have a 
degree less or equal to 6, and therefore the first stage of the 7- Color-Planar- Graph 
algorithm terminates in at most O(log n) steps. At each step we have to identify the 
vertices that have degree less than 7 in the remaining graph. This takes constant 
time on a CRCW PRAM (assuming that if two or more processors simultaneously 
write into some location, one of them will succeed) and O(log n) time on an EREW 
PRAM. 

In the second stage all the uncolored vertices are of degree less or equal to 6 and 
therefore, by theorem 4.4.2, the procedure Constant- Degree- MIS finds, in 0(lg*n) 
time, an MIS in the graph induced by these vertices. When the algorithm colors a 
maximal independent set, at least one uncolored neighbor of each uncolored vertex 
becomes colored. Therefore the second part of the second stage terminates in at 
most 7 iterations. 

Since the first stage takes 0(log n) time on a CRCW PRAM and 0(log 2 n) time 
on an EREW PRAM, and since each one of the O(log n) iterations of the second 
stage is dominated by a call to Constant- Degree- MIS ', the total running time is 
0(log n lg*n) on a CRCW PRAM and O(log 2 n) on an EREW PRAM. I 

Remark: If, at each stage, instead of removing from the graph all the vertices with 
degree less than 6, we remove all the vertices with degree less or equal to the average 
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degree, the algorithm described above produces a correct result in polylogarithmic 
time for any graph G such that the average degree of any vertex-induced subgraph 
G' of G is polylogarithmic in the size of G'. This class contains many important 
subclasses including graphs that are unions of a polylogarithmic number of planar 
graphs (i.e. graphs with polylogarithmic thickness). 

Our techniques together with the ideas presented in [10] can be used to construct 
a deterministic O(log 3 n lg*n) time algorithm for 5-coloring a planar graph. The 5- 
coloring algorithm has two stages. The first stage of the algorithm partitions the 
graph into layers such that vertices in any layer are independent and have degree of 
at most 6 in the graph induced by the vertices in its layer and the higher numbered 
layers. This partitioning can be done by partitioning the graph into layers as in the 
first stage of the 7-coloring algorithm, coloring each layer in 7 colors, and taking 
each color class to be the new layer. The second stage of the algorithm adds layers 
one by one, starting from the layer with the highest number, each time recoloring 
the graph with 5 colors. 

Before describing the second stage, we need the following definitions. Let G be a 
partially colored graph and let c x and c 3 be two distinct colors. A color component 
is a connected component of a subgraph of G induced by all vertices of color c\ and 
C2. A color component flip is a recoloring of the color component that exchanges 
colors c\ and c 2 . A color component flip does not affect the validity of coloring. 

We can proceed with the description of the second stage of the algorithm. After a 
layer is added to already colored graph, we first color all vertices that can be colored 
without changing the existing coloring. This can be done in the same way as in 
the 7-coloring algorithm. Now all 5 colors are represented among neighbors of each 
uncolored vertex. Since the uncolored vertices have degree of at most 6, the results 
of [10] imply that for every uncolored vertex v there are two colors c\ and cj such 
that t; has exactly one neighbor w\ of color c\ and exactly one neighbor u> 2 of color c 2 . 
Furthermore, the vertices wi and u>2 belong to different color components induced 
by colors C\ and c<i. Flipping each one of these color component allows us to color v. 
The problem is, however, that flipping both color components simultaneously does 
not allow us to color v. We call such color components dependent. 
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Where as Boyar and Karloff use randomness to deal with these dependencies, we 
use our symmetry-breaking techniques as follows. For each pair of distinct colors c x 
and C2, we construct color components induced by these colors. Then we construct 
a dependency graph with vertices corresponding to the color components and edges 
corresponding to the dependencies between the color components. Flipping a set of 
color components that corresponds to an independent set in the dependency graph 
does riot cause conflicts. Suppose we can find an independent set in the dependency 
graph such that nipping the corresponding set of color components allows us to color 
a constant fraction of uncolored vertices. Then in O(log n) iterations will be able 
to color all uncolored vertices. 

We find such an independent set in the dependency graph as follows. Observe 
that the dependency graph is planar, so we can 7-color this graph using the 7- 
Color- Planar- Graph algorithm. Then, for each pair of distinct colors and for each 
color class of the corresponding dependency graph, we compute the number of 
uncolored vertices of the original graph which can be colored if the color components 
corresponding to vertices in the color class are flipped. For each of the 10 possible 
choices of colors c\ and c 2 there are 7 color classes, so the total number of times 
that we count the number of vertices that can be colored if a color class is flipped is 
70. Since each uncolored vertex is counted at least once, there is a color class such 
that flipping all color components in this class allows us to color at least 1/70 of all 
uncolored vertices. 

Now we analyze to complexity of the algorithm. The outer loop of the algo- 
rithm that iterates over layers is executed 0(log n) times, and the inner loop that 
colors a constant fraction of uncolored vertices is executed 0(log n) times as well. 
Each iteration of the inner loop does 10 connected component computations, 70 
enumerations, and 10 calls to the 7- Color- Planar- Graph procedure. Since a con- 
nected component computation can be done in O(logn) time on CRCW PRAM 
using Shiloach-Vishkin algorithm [65] and an enumeration can be done in 0(log n) 
time by using parallel prefix computations (see Section 3.2), the 7- Color- Planar- 
Graph procedure, which runs in 0(log n lg*n) time is the bottleneck of the inner 
loop. The overall running time of the algorithm is 0(log 3 n lg*n). 



4.5. ALGORITHMS FOR PLANAR GRAPHS 111 

The above result is summarized in the following theorem. 

Theorem 4.5.2 A planar graph can be 5-eolortd in 0(k>g 3 n) time on a CRCW 
PRAM using 0(n) processors. 

Using the techniques described in this chapter it is easy to construct a fast 
algorithm for finding a maximal matching in a planar graph. 

Theorem 4.5.3 A maximal matching in planar graph can be found in 0(log n lg*n) 
time on a CRCW PRAM. 

Proof: First, the algorithm partitions the graph into 0(lg n) layers, such that the 
vertices in each layer are independent. This is done by iteratively removing a 
graph induced by vertices of degree six or less, 7-coloring the removed graph, and 
considering each color class to be a layer. The algorithm proceeds by iteratively 
adding a layer, finding a maximal matching in the obtained graph, and removing 
the end-points of the edges in the matching. 

At the end of each iteration the remaining set of vertices V\ induces a graph 
of degree zero. Next a new layer with the vertex set Vj is added. By construction 
of the set Vj, the degree of each added vertex is at most 6. Therefore we get a 
bipartite graph with vertices in Vi on one side and vertices in Vj on the other side. 
Maximal matching in this graph can be found as follows. Each vertex in V 3 marks 
an adjacent edge. Each vertex in Vi chooses one of its marked adjecent edges. The 
chosen edges are added to the matching and their end-vertices are removed. It is 
clear that after 6 iterations all edges are removed. 

Each iteration over a layer takes 0(lg*n) time on a CRCW PRAM and the 
number of iterations is 0(log n). This gives 0(log n lg*n) total running time. I 

Proof: First, the algorithm partitions the graph into layers, such that the vertices 
in a layer are of degree at most 6 in the graph induced by the vertices of this layer 
and higher-numbered layers. Each layer has degree of at most 6 and therefore can 
be 7-colored in 0(lg*n) time. The algorithm proceeds by iteratively adding a layer, 
finding a maximal matching in the obtained graph, and removing the end-points 
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of the edges in the matching. At the end of each iteration the remaining vertices 
induce a graph of degree zero, and these vertices are independent in the graph 
existing at the current iteration. Since the graph added at the beginning of the 
current iteration is 7-colored, by coloring the remaining vertices with a new color 
we obtain an 8-coloring of the current graph. This coloring gives us a coloring of 
the edge-graph of the graph using 56 colors. Hence, a maximal matching in the 
current graph can be found in constant time by finding a maximal independent set 
in the edge-graph. Each iteration takes 0(lg*n) time on a CRCW PRAM and the 
number of iterations is O(logn). This gives O(lognlg*n) total running time. | 



4.6 Lower Bounds 

In this section we prove two lower bounds for a CRCW PRAM with polynomial 
number of processors: 

• Finding a MIS in a general graph takes ft(log n/ loglog n) time. 

• 2-coloring a directed list takes J2(log n/ loglog n) time. 

The first lower bound complements the O(log n) CRCW PRAM upper bound 
for the MIS problem that is achieved by Luby's algorithm [55]. The second lower 
bound complements Theorem 4.3.2 in this chapter. 

We need to define the following two functions. The PARITY function on n bits 
xx . . . x„ equals {x x + . . . + x„)mod2. The MAJORITY function on n bits equals 1 
if xi + . . . + x n > n/2, and equals to otherwise. The results of [25,4,5] imply that 
any CRCW PRAM algorithm for PARITY or MAJORITY that uses a polynomial 
number of processors requires Q(log n/ log log n) time. 

Theorem 4.6.1 The running time of any MIS algorithm on a CRCW PRAM with 
a polynomial number of processors is fi(logn/ loglog n). 



4.6. LOWER BOUNDS 113 

Proof: Given an instance of MAJORITY, we construct an instance of MIS in con- 
stant CRCW PRAM time. Let x 1 ,x 3 ,...,x n be an instance of MAJORITY. We 
construct a complete bipartite graph G = ( V, E) with vertices corresponding to zero 
bits of the input on one side and vertices corresponding to one bits on the other 
side. The graph G = ( V, E) is defined by 

V = {l,...,n} 

E = {(i,j)\xi^xj}. 

To construct this graph, assign a processor P,, for each pair 1 < * < j < n. Then, 
each processor Py writes 1 into location M y if x, / Xj and otherwise. 

A maximal matching in a complete bipartite graph is also a maximum one. By 
constructing a maximal independent set in the line-graph G' of G, one can find a 
maximal matching in G. To construct the graph G' assign a processor P tJ * for each 
distinct i,j, k < n. Each P ijk writes 1 into location Af(<j),(j,fc) if My = M ik = 1 and 
otherwise. 

The MAJORITY function equals 1 if and only if there is an unmatched vertex 
* € G such that Xj = 1, which can be checked on a CRCW PRAM in constant time. 
I 

Theorem 4.6.2 The time to 2-color a directed list on a CRCW PRAM with a 
polynomial number of processors is fi(logn/loglogn). 

Proof: We show a constant-time reduction from PARITY to the 2-coloring of a 
directed list. First, we show how to construct, in constant time, a directed list with 
elements corresponding to all the input bits x< with value of 1. Let xi, X2, . . • , x„ 
be an instance of PARITY. Associate a processor Pi with each input cell Af,- that 
initially holds the value of x,-. Associate a set of processors P/ with each index 
*> 1 <• k < j < i. In one step, each processor P/ reads the value of Mk and, if 
it equals 1, writes 1 into M/, effectively computing the OR-function on the input 
values x<_j, x,_j + i, . . . , x t _i. Assign a processor P/ to each Mj . Each processor P/ 
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reads Mf and M/ +1 and writes j into M/ if and only if M/ ^ M/ +1 . It can be seen 
that for all < t < n, M[ holds max{j | ,;' < i,Xj = 1}. 

We have constructed a directed list with elements corresponding to all the input 
bits Xi with value of 1. Assume this list is 2-colored. Then PARITY equals 1 if and 
only if both ends of the list are colored in the same color, which can be checked in 
constant time. | 



4.7 Remarks 

In this chapter we described the technique for 3-coloring rooted trees and showed 
how this technique can be used to design parallel algorithms. Our technique applies 
when we have a set of jobs each of which can be performed at this time, but 
simultaneous execution of these jobs may cause conflicts. As we have seen, if the 
graph describing these conflicts is sparse, the technique allows to select a large 
subset of non-conflicting jobs. Note that there is no need for such a technique in 
the sequential framework. This chapter shows how important it is to identify the 
fundamental problems of parallel computation and to developed efficient solutions 
to these problems. 



Conclusion 

In the first chapter of this thesis, we introduced a new approach to the maximum 
flow problem. This approach combines the concept of preflow due to Karzanov 
with a novel concept of distance labeling. Using our approach, we have developed 
better algorithms for the problem, including an 0{nm log(n 3 /m)) time sequential 
algorithm and 0(n 2 log n) time parallel algorithm. The parallel algorithm uses a 
linear amount of memory and, as we have seen in Chapter 3, is practical. 

In the second chapter, we developed a framework for solving minimum-cost 
flow problems. This framework permits the extentions of maximum flow tech- 
niques to the more general minimum-cost flow problem. We generalized the tech- 
niques developed in Chapter 1 to obtain an 0(nm log(n) k>g(nC)) time sequential 
algorithm. We also believe that a generalization of these techniques leads to an 
0(mnlog(n 2 /m)log(nC)) time sequential algorithm and an 0(n 2 log(n) log(nC)) 
time parallel algorithm that uses a linear amount of memory, but the bounds we 
prove in Chapter 2 are somewhat weaker. We also show how to incorporate in our 
framework Dinic's blocking flow method for solving maximum flow problems. Using 
the algorithms for finding blocking flows in acyclic networks, we obtain sequential 
algorithms for the minimum-cost flow problem that run in 0(nmlog(n)log(nC)), 
0(n 5 / 3 m 2 / 3 log(nC)), and 0(n 3 log(nC)) time, and an 0(n 2 log(n)log(nC)) time 
parallel algorithm that uses 0(n 2 ) memory. 

In the third chapter we described an implementation of the parallel maximum 
flow algorithm. The key feature of this implementation is the use of parallel prefix 
operations as primitives, which enables us to obtain an efficient implementation 
while maintaining a high level of abstraction. The usefulness of the parallel prefix 
operations stimulated our search for other fundamental techniques for the design 
and implementation of parallel algorithms. In Chapter 4 we introduced parallel 
symmetry-breaking techniques and showed how these techniques can be applied to 
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design efficient algorithms. We also showed lower bounds which imply that our 
symmetry-breaking techniques are optimal in some contexts. 

In conclusion, we would like to discuss potential improvements and extensions 
of the results of this thesis. As discussed in Section 1.7, an interesting question 
related to the maximum flow problem is the existence of an 0(nm) algorithm for the 
problem. For the special case of zero-one networks, Even and Tarjan [17] show that 
Dinic's algorithm is more efficient than for the general networks. Can their bounds 
be improved using the approach of Chapter 1? The improved bounds would imply a 
better bound for the maximum bipartite matching problem. For the minimum-cost 
flow problem, a promising approach is scaling by value (i.e., cost-capacity product), 
as discussed in Section 2.11. 

The maximum flow problem and the minimum-cost flow problems have been 
studied extensively. The multicommodity flow problem, on the other hand, has 
received relatively little attention. The current algorithms for this problems [43,45] 
use the techniques of mathematical programming. An efficient multicommodity flow 
algorithm that uses graph-theoretic techniques similar to the techniques discussed 
in this thesis would be a very interesting result. 

In the area of parallel algorithms, one of the most important problems is that of 
identifying the fundamental problems related to parallel computation and finding 
efficient solutions for these problems. 

This thesis makes a step towards a better understanding of algorithms on graphs, 
both in sequential and parallel contexts. It is our hope that the techniques and 
approaches presented here will prove useful in further advances in this important 
field. 
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W - design efficient algovtfhffle for the-miaimaawcaat flow PXQbJ.em* „,«,Mff, exhibit 
< ' 0(nmlog(n)log(nC)), OCnS^^^Kirf)^*^;^^^*^ time sequential 

|r* algorithms as well as ya&glxk v*4 &s*&b**j*^+lw^tim*^ -• : . 

I ' In Chaptef 3 we address inp$«neirfcatlon of- parallel -algorithms through a 
f -case-study of an impleite&tation of a parallel maximum flow algorithm. 

Parallel prefix operations play an, important "TOle in crff ^plementsticm. 
I We present experl^ntal^results achieved! by the implementation. 

Parallel symmetry-breaking techniques are the main topic of Chapter 4. We 
I give an 0(l £ *n) algorithm for 3-coloringi a rooted tree. This algorithm is 

! "' ' 'used' to improvte e##»^a^ilAel\*aa^itJSn&^ »>*,**_ 1- "■< 

! coloring, and finding, maSdmal independent^ set in coiistanf-degree graphs, 5- 

l coloring Planar graphs, and finding a makimal matching In planar graphs. We 

also prove Tower bZZZU oSI t^ >raliet g-ompTekif y oT 'fnT ffilbei^t . independent 

set problem and the problem of 2-coloring a rooted tree. 
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